Modelling of crust and mantle heterogeneity effects in the theory of Earth tides Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 6, NO. 2, PAGES 80-93, APRIL 2004

Modelling of crust and mantle heterogeneity effects in the theory of Earth tides

S. M. Molodensky

Schmidt United Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia

Abstract. General relations are obtained for the determination of radial and tangential components of displacements and gravitational potential variations at the Earth’s surface due to tidal forces for a model with a horizontally inhomogeneous distributions of elastic moduli. Model calculations are performed to estimate the effects of local and regional inhomogeneities of elastic moduli specified in regions of simple configurations on tidal variations in gravity, tilts and deformations. As distinct from gravity data, data on tidal tilts and deformations can provide significant constraints on local inhomogeneities of elastic moduli in the immediate vicinity of an observation point. Results of model calculations show that the phases and amplitudes of tidal tilts and deformations change by 10-15% in the vicinity of a soft inclusion in which the shear modulus is smaller than its value in the surrounding medium by 20-30%. Anomalies of such a contrast can be detected over several months of continuous tidal observations.

Introduction

The study of local heterogeneities in the crust and upper mantle by traditional seismic methods is largely complicated due to uncertainties in the parameters of earthquake sources. The method of controlled vibratory sounding is free from this problem, but its application requires expensive sources of seismic signals, and their effective application is often inappropriate. The lunisolar tides are excited by natural sources of tidal forces, whose values are known with a reasonably good accuracy at all time moments. Therefore, in many cases the use of data on the material response to tidal forces seems quite advantageous for the study of mechanical properties of the medium. We should also note that tidal data differ from seismic evidence in that they provide information on the properties of the medium in the range of much lower frequencies (diurnal and semidiurnal periods). Because of this distinction, comparison of seismic and tidal data can provide constraints on both the elastic moduli and

Copyright 2004 by the Russian Journal of Earth Sciences.

Paper number TJE04149.

ISSN: 1681-1208 (online)

The online version of this paper was published 3 June 2004.

URL: http://rjes.wdcb.ru/v06/tje04149/tje04149.htm

rheological properties in the vicinity of an observation point.

In recent years, the most significant anomalies of tidal tilt and strain factors have been discovered in the vicinity of the Orenburg gas-condensate field and in the Baksan Gorge at the foot of the Elbrus Mountain. In both cases, this appears to be related to the presence of a cavity (possibly containing a rigid skeleton) significantly differing in average shear modulus from the surrounding medium. Analysis of the response of this cavity to tidal forces allows one to estimate (1) the sizes of this cavity and mechanical properties of its upper dome and, based on the estimated thickness of the dome and its effective rigidity, (2) the criteria of its stability (this problem is very important in relation to the possible collapse of the dome due to the pumping-out of gas).

Along with the detection of local heterogeneities, there also exists the problem of identifying global lateral heterogeneities in the mantle, which should primarily give rise to anomalies in the tidal gravity factor S (its value is mainly dependent on mechanical properties of the lower mantle and nearly insensitive to local heterogeneities in the vicinity of an observation point).

Below, we derive general relations for the determination of radial and tangential components of displacements and variations in the gravity potential at the Earth’s surface with laterally inhomogeneous distributions of elastic moduli. Numerical estimates are presented for the effect of local inhomogeneities of elastic moduli on tidal tilts and deformations.

1. Initial Equations

Partial differential equations in a general form governing the elastic-gravitational conditions of the Earth under the action of tidal forces are based on the assumption of an initial (nondeformed) hydrostatic equilibrium of the mantle and crust, implying that the pressure gradient Vp0 is compensated for by gravity:

Li(u, R)

(henceforward, summation is assumed under repeated indexes and no distinction is drawn between covariant and contravariant components of vectors and tensors, i.e. all calculations relate to a Cartesian system of coordinates).

Vpo - pVW = 0 , (1)

where p is density; W is the unperturbed gravitational potential, satisfying the Poisson equation

AW = -4nGp , (2)

and G is the gravitational constant.

Elastic tidal deformations displace an element of the medium from a point determined by the radius vector r to the point r+u and compresses this element, with the relative St

change in volume being — = V • u. Then, the density p',

T

potential W' and stress tensor a'k in the deformed medium can be expressed as

p'(r) = P(r - u) - PV • u = P(r) - V • (Pu) , (3)

2. Generalization of the Betti Theorem of Reciprocity to the Case of Tidal Deformations of a Medium.

Green Function

To describe tidal deformations in a radially and laterally heterogeneous medium, we use Equations (8) and (7b) with the inhomogeneous boundary conditions

a • n| =0 , (9a)

W'(r) = W(r) + R(r) ,

Sik (r) = -po(r)Sik + Sik (r)

= [-po(r) + (u, Vp)]Sik + Sik

(4)

(5)

_ ^ , (dui , duk\

sik = A(V • u)Sik + № ydx + J

((OR/ dn) int1 s - (dR/dn)out) = 5(vt(a)/a)y2('&,$,t) ,

(9b)

where Sik is the Kronecker delta and sik is the tensor of elastic stresses connected with the displacement vector through the Hooke law

(6)

Substituting (3-6) into the equation of elastic equilibrium and subtracting, respectively, (1) and (2) from the results, we obtain

pVR - (V • pu)VW + V(p(u, VW))=0 . (7a)

where E and n are the Earth’s surface and its outer normal; vt is the amplitude of the tide-generating potential; a is the average radius of the Earth; y2(’&, f, t) is the spherical function of the 2nd order; '&, f and t are colatitude, longitude and time, respectively; and the indexes “int” and “out” indicate, respectively, inner and outer normal derivatives.

To calculate the Green function for the boundary problem (8), (7b), (9a), (9b), we will write relations generalizing the Betti theorem of reciprocity to the case of tidal deformations of a gravitating medium with a hydrostatic distribution of initial stresses.

In order to solve an inhomogeneous equation of the form

L(u(1),R(1)) = F(r)

(10a)

AR = 4nGV • (pu) .

(7b)

AR(1) = 4nGV • (pu(1))

(10b)

where F(r) is a known function of spatial coordinates, we

Applying the curl operator to the left- and right-hand sides

introduce the following system of auxiliary solutions (u(j), R(j)) of homogeneous equations (8) and (7b):

L(u(j),R(j))=0 , (10c)

of (1) yields

VW x Vp = 0 .

Adding the identically zero term

u x (VW x Vp) = (Vp, u)VW - (VW, u)Vp ,

to the left-hand part of Equation (7a), this equation takes the form

AR(j) = 4nGV • (pu(j)) which satisfy the boundary conditions

(10d)

s

„(j)

n|s= C1erym(^,0) ;

D

(j)

0 at j = 1 ;

(11a)

solution (u(j), V(j))); the additional terms

J [D1R(j) - D(j)R1]da

s(j)

D(j) = C2ym(^,0) at j = 2 ;

0

„(j)

D

s(j)

(j)

|s= C3erVym(tf,0) ; 0 at j = 3 ;

n|s= c4er x Vym(tf,0) ;

D(j) = 0 at j = 4 .

(11b)

(11c)

(11d)

Here

I = j [((S1 • n), u(j)) ((s(j) • n), u1)

s

+ D1R(j) - D(j)R1]da ,

determine the difference between the gravitational energies of the single layer D1 in the field with the potential R(j) and the single layer D(j) in the field with the potential R1.

If the external body forces F(r) are nonzero, relations (13a) and (13b) yield the identity

(15)

D(j) = ((9R(j)/dn)int|s-(dR(j)/dn)out |s)/(4nG) (12)

is the density of the single layer on the Earth’s surface, er is the unit radius vector and yjf are normalized spherical functions.

Taking into account the explicit form of the operator L (8), it is easy to show that this operator is Hermitean, i.e. the volume integral of scalar product

1 = /(^^ L(ul,Rl))-(u(1), L(u(j),R(j))))dT (13a)

reduces to the surface integral

I = y (u(j), L(u1,R1^dT

r

= j (u(j), F(r))dT = j [((s1 • n), u(j))

t S

- ((s(j) • n), u1)+D1R(j) - D(j)R1]da ,

which is helpful for the construction of the Green function related to the problem of elastic deformations of a gravitating medium produced by the given external body forces F(r).

Indeed, if external surface forces are absent, we have s1 • n = 0. Let the solutions (u(j), V(j)) be specified by the boundary conditions s(j) • n = e(j)S(r - r0) and D(j) = 0, where e(j) is an arbitrarily oriented vector, S(r - r0) is the 2-D delta function, and r and r0 are vectors on the Earth’s surface. Then, integrating the right-hand part of (15) over angular variables, we obtain the relation

(u1 (ro), e(j))=y (u(j), F(r))dT

(16)

(13b)

where the density of the single layer D1 is determined by a relation similar to (12):

D1 = ((dR1/dn)int|S-(dR1/dn)out|S)/(4nG) . (14)

If we set F(r) = 0 in (10a), i.e. if the solution of homogeneous tidal equations in the absence of external forces is taken as (u(1),V(1)), volume integral (13a) vanishes identically, and relation (13b) generalizes the Betti theorem of reciprocity to the case of a gravitating medium with a hydrostatic distribution of initial stresses. The terms

j [((s1 • n), u(j)) — ((s(j) • n), u1 )]da s

in (13b), corresponding to the theorem of reciprocity, determines the difference between the values of the work done by surface forces acting through the displacements u(j) (corresponding to the solution (u(1), V(1))) and by surface forces acting through the displacements u(1) (corresponding to the

which completely defines the three components of the sought-for vector u1 at any point ro of the surface a.

Analogically, specifying the solutions (u(j),V(j)) by the boundary conditions s(j) • n = 0 and D(j) = S(r - ro) and integrating the right-hand part of (15) over angular variables, we obtain the relation

R(j)(ro) = j(u(j), F(r))dT , (17)

which completely defines the potential variations on the Earth’s surface.

3. Application of the Perturbation Method to the Problem of the Effect of Horizontal Heterogeneities on Elastic Deformation of a Medium

Tidal displacements and the potential in spherically symmetric models of the Earth are determined by the Love formula

n

s

n

uo = f1o)(r)eryk + '~6 f2o)(r)Vyk , (18)

Ro = /3o)(r)y2k ,

where k = 0,1 and 2 for long-period, diurnal and semidiurnal tidal components, respectively.

When a horizontally homogeneous model of a medium (with distributions of density po(r) and elastic moduli Ao(r) and ^o(r)) is replaced by a spherically nonsymmetrical model (with distributions po(r) + Sp, Ao(r) + SA and ^o(r) + S^), the solutions (uo,Ro) of Equations (8) and (7b) acquire increments (u1,R1).

In the first approximation of the perturbation theory, the equations of tidal strain in a laterally and radially heterogeneous medium can be written as

u(j) = t(j)(r)er xVym at j =4 R(j) = 0 .

(21b)

Here fi (r) is a system of six functions defined by the ordinary differential equations

Lo(u1, R1) + L1(uo, Ro) = 0

AR(1) = 4nGV • (pou(1) + p1u(o))

where

L1 (uo, Ro)

Sp(V(Ro + (uo, VW)) - (V • uo)VW)+V • s1

and

u(1) = E

S/1nm)(r)er ym +

t(n + 1)

T _f(n,m) / \ v"7

S/2 )(r)Vy

/i = E aik /k

(19a)

(19b)

(19c)

with the coefficients

an = -a66

a12 = -a56 =

a16 =

2A

r(A + 2^) A ;

A + 2^

r2(A + №)

-a21 = a65

-n(n + 1)

a31 = -a64 = 4nGp ;

n + 1

-a33 = a44 = --------- ;

a34

4nG

s1ik = SA(V • uo)Sik + S^ddu^ + . (19d)

The general solution to the problem of tidal strain in a horizontally heterogeneous medium is sought as the sum of spheroidal and toroidal spherical harmonics:

a52 = 2^r'

f 2(A + ^) 1

a61 = 4prW + 4^

y A + 2^ n(n + 1)

3A + 2^

A + 2^

(22)

(23)

a41 = a63 = (n + 1)pr ;

2

a42 = a53 = -pr ;

2™ 7' r, 3A + 2№

a51 = a62 = -pr W - 2^r ;

A + 2^

+ ]T [St(nm)(r)er x Vym]

(20a)

The symmetry properties of the coefficients aik in (23) are due to the self-conjugacy of the operator L (8) [Molodensky, 1984].

The boundary conditions for the solutions (u(j), R(j)) are defined as follows:

r(1) = E S/3nm)(r)ym.

(20b)

Auxiliary solutions (u, R(j)) determining the Green function are represented as expressions of the same form as the terms in (20a) and (20b) at j = 4 :

2

u(j) = /1j)(r)erym+ „J, /V^ym at j = 4

n(n + 1)

R(j) = /3(j)(r)ym, at j = 1, 2, 3;

and

(s(j), n) = erynm; D(j) =0 at j = 1; (s(j), n) = Vynm; D(j) =0 at j = 2; (s(j), n)= 0; D(j) = ynm at j = 3;

I = S/1nm) at j = 1 ;

I = n(n + 1)S/2(nm) at j = 2 ;

I = S/3nm) at j = 3;

I = t(nm) at j = 4 .

(21a)

Here

(24)

(s(j), n) = er x Vynm; D(j) =0 at j = 4 . Then, substituting (20) and (21) into (13), we find

1

2

r

r

2

r

2

r

I = / (u(j), L1(uo,Ro))dT .

(26)

I=

(u(j)) ^) dT :

v 7i dxk )

(du(j))i i d

s1ik—^------- I dT

dxk

These relations completely define all of the sought-for coefficients in expansions (20a) and (20b).

4. Green Function for the Problem of Tidal Deformations in a Heterogeneous Medium

In order to numerically calculate the Green function in a heterogeneous medium, we express the components of the tidal displacement vector through homogeneous harmonic polynomials £nm = rnyjm($, f). In Cartesian coordinates xi, they satisfy the equation

The substitution of (19d) into (20) yields I = I(SA) + I (S^) ,

where

(31a)

and

I (SA) = -J SA(V • uo)(V • u(j) )dT (31b)

7VT 1 /\ ((duo)i , d(uo)^ (du(j))i

(S№) = - I SM + -j— ) -axr-dT . (31c)

d£n

dxi

= n£„

(27) Substituting (29a) and (29b) into (31c), we obtain

and the Laplace equation

((duo)i + d(uo)^ (du(j))

A£nm = 0 .

(28)

y dxk dxi J dxk

= h2£o1 {2r3hnj) £nm + 2r2hnj)Cn r

Expressing the displacement vectors uo and u(j) through £nm, we obtain

(j)

+ 2r2n£nm I hn ) + — ) + 2tnj)Xi

d2£ 1

nm

(uo)i = h(o)xi^o + t(o) ^

(n(j)V = fr(j)r-£ + t(j) —

(U )i — hn xiSnm + t-n o

dxi

(29a)

(29b)

where £o = r yj($,f); j = 0, 1 and 2 for long-period, diurnal and semidiurnal waves, respectively; and the functions hnj),tnj) and h(o),t(o) depend on the radius r alone. Comparing (18), (21a) and (29), these functions can easily be expressed through /i(j) (r), determined by the system of equations (22)-(23).

Now, we address the solutions of tidal equations in the case of arbitrary lateral variations in elastic moduli. Then, integral (26) takes the form

n xixk n n f

dXidXk J

+ fo£o j 2rhnj) £nm + +6hnj)Cnm

+ 2n£nm ^hnj) + |

, (, , ^2 ) ( d£o , d£o \

+ h2 +-------( xit;-------+ xk 7— )

y r J V dxk dxi/

x { hj) £nm + hnj)^nmSik

+ 1 ('*' + j )(x. if + Xk-

(32)

d^nn

dxi

+t + t2

(j) d2£n

dxidxk ^

£o lnL(j)' XiXk £ + 2h(j)£ S't

Snm + 2hn £nmSik

d2£o

dxidxk

I=

(j)

ds1

dT

where

_ X\^V7 ^ , X ( d(uo)i d(uo)k

s1ik — SA(V • u)Sik + S^ ( ——--------------+

dxk

dxi

(29c)

(29d)

+

(h?> + j (

. tnjW d£nm , d£nm

+------*i^------+ Xk 7:

r V dXk dxi

(j) d2£nm

+ tn dxidxk Taking into account (27), we have

Using the Gauss theorem and taking into account the condition of vanishing stresses at the outer surface,

dxk dxidxk dxi V dxik dxk

s1iknk = 0

we obtain

d£n

(xk A£o + dxoSi^ = (n - 1)(V£nm, V£o)

dxk

2

nm

t d2£nm _ d f d£o c d£o \

xixk£nm^ ^ — TJ ( TJ xixk£nm TJ )

dxidxk dxi Vdxk dxk/

nm

x ^xixk _d — + 3xk£nm + xiSik£n

(33)

d

— d (2£oxi£nm) 2 (n + 4) £o£nm — 2£o£n

a2(r) = (r)r l j ^hnj) + tir_^ ^ho +

+ 2(n - 1)tnj^ho + + 2to ^hnj) +

a3(r) = 2№;P(r)r-itn7)to .

The expansion of the lateral inhomogeneities of elastic mod- For calculating the integrals of products of three harmonic uli in series of harmonic functions yields polynomials, we substitute the relation £nm = rnym ($,f)

into the first term in (33b) and obtain

SA(r,^,f) Aip(r)y p(^,f) ,

l = 1 p=-l oo l

SM(r,^,f) = EE №ip(r)y p(^,f) .

(33a)

£ip£nm£oa1(r) j dT

(34)

l = 1 p=-l

Ai

pnmngmg

r

l+n+no + 2

a1(r)dr .

Taking into account (32)and (33), the integrals in (31b) and o

(31c) can be represented as the sum of terms containing in- Here Aipnmncmc is the integral of the product of three spher-

tegrals of products of three harmonic polynomials, scalar ical functions with the respective indexes (1,p), (n,m) and

products of their gra^^ and the convolution of the ten- (no, mo) taken over the surface of a unit sphere.

sors consisting their second derivatives and having the form The second term in formula (33b) contain the product

of the homogeneous harmonic polynomial £lp and the scalar product of the polynomial gradients (V£nm, V£o) and can be calculated as follows. Denote

I(SA) = - J SA(V • uo)(V • u(J))dT

T

= E JJ j {£lp£nm£o (r)C(r)|dT ,

l,p t

-I(S№)= f Sfj, ((duo)i ■ d(uo)k A (du(j))

(1)

nomonm

£lp( V£nm , V£nomJ dT . (35)

+

dxi

V dxk

T

^ ' III \ £lp£nm£oa1(r)+ £lp (V£nm, V£o) a2 (r)

dT (33b) Integrating (35) by parts and taking (27) into account, we

find

J(1) + J n0m0 = „a

Jn0m0nm + Jlpnm nAIpnmnomo

(36)

+ £lp

d2£nm d2£o

dxidxk dxidxk

«3 (r) > dT

and two similar relations obtained from (36) by the cyclic permutation of the polynomials £lp, £nm, £nomo. Solving these equalities, we find

where

Jnomonm 2 (n + no ^) Ai

ipnmnomo

(37)

Z(r) = Aip(r)r-(l+n+2) ((/(j))'(r) + 2/(j)(r)

- f2(j) (r))((/l(o))'(r) + 2 A(o)(r) - f2o)(r)) ,

a 1 (r) = 2^p(r)r-i jr2hnj) h'o + r(hj) ho + h'ohn^) (j)'

The third term in (33b), containing the volume integral of the convolution of tensors that consist of the second derivatives of homogeneous harmonic polynomials of the form

£ij

d2£nm d2 £o

dxidxk dxidxk

+ nrh'o ( hnj) + —— I + 2rhnj) ( ho + —

(33c)

can be calculated in a similar way. Denoting

+

n(nr 1) h/otnj) + 2hnj) to + 3hnj)ho

+ nho | hnj) + — | + 2hnj) fho + — rr

jlp(2)

no mo nm

£ij

d2 £nm d2 £o

dxidxk dxidxk

dT

(38)

and integrating (38) by parts, we obtain

jlp(2) __ /*/* £ d£nm d £o , jnomo(3)

Jnomonm = xi£lp "^7 dxidxk ds - Jlpnm , (39)

where

2

r

R

r

r

Jn0m0

ipnm

(3)

d£ip d£nm d2£o

dxi dxk dxidxk

dT

+ 4 [(l + n + no)2-1](n + no - l)(n + no - l - 2) ar2r) j

(40)

Taking (27) into account, we have

d2£o = / d£o\ _ d£o

dxidxk dxk V dx^ dx

_ ( -1) d£o

= (no - ^ dik

Sik

(41)

and, therefore,

*i£l;

d£nm d2£o

p dxk dxidxk

ds

(no - 1H £lp(V£ nm, V£o )ds

Jlp(1) nmn0m0

= (no - 1) (l + n + no + 1) Jn

= 2 (no - 1)(n + no - 1)(1 + n + no + 1) Alpnmn0m0

Integrating J-:mo(3) by parts, we obtain

rnomo(3K Tnm(3)

J + Jl

lpnm

d£lp d£nm d£o ,

xi ^ ----t:—ds

1 + n + no + 1) J,

lp(1) ( nomonm (

no mo lpn

(3)

1

2 (1 + n + no + 1)

nomo(1) + i jlp(1)

( — „n Jn0m0(1) + J Jlp(1) + nJnm(1) \

V noJlpnm + 1Jnmnomo + /Jlpnomo/

or, substituting (37),

Jrp0T(3) = 1 (1 + n + no + 1)

x („o - n2 - 12 + 2nZ) Ai pnmnomo

The substitution of (42) and (45) into (39) yields

Jiip(^)[inm = 1 (l + „ + no + 1)(no + n - l)

(no + n - 1 - 2) Ai

pnmn0m0

With the help of (37) and (36), expression (32) is transformed into the final form

-I(S^) = E Alpnmn0m0 J drrl+n + n0 ^ (r2 «1 (r) l,p o

+ 2 (l + n + no + 1)(no + n - l) «2 (r)

5. Tidal Deformations of a Heterogeneous Elastic Half-Space

In analyzing local effects of heterogeneity of a medium on tidal tilts and strains, it is sufficient to consider the simplest model of an elastic and nongravitating half-space.

We consider a system of two auxiliary solutions (corresponding to the values j=1 and 2) of the homogeneous equations

Li (■

(j) (j)

V(j) = 0, at j = 1

(42)

AV(j) = 4nGV • (pu(j)) at j = 2

(47)

(48)

with boundary conditions of the delta function type at the outer surface (corresponding to the cases of tangential and normal concentrated loads applied at a point (xo,yo, 0):

lpnomo i i "i dxi g^k (43)

Applying the cyclic permutation of indexes (1,p) ^ (n,m) ^ (no, mo) to this relation and determining JlpOf from the resulting system of three equations, we find

a(j)nk = aoS(x - xo)S(y - yo)S(z)ex at j = 1 ,

aj)nk = aoS(x - xo)S(y - yo)S(z)ez at j = 2 (

(49a)

(49b)

(44)

Here (x,y,z) is a local Cartesian system of coordinates oriented in such a way that the z axis is directed along the normal to a surface element and the x axis coincides in direction with the action of load in the j = 1 solution.

Substituting (47) and (48) into (49a) and (49b), we obtain

(45)

aou“,(w’’0) — /<u“

T

= - J SA[(V • u(j))(V • uo)+2S№elk)eiko]dт ( )

(j = 1)

(46)

(1)

aouZ

(xo, 0)= J (u(j), F(r))dT

T

= - J SA[(V • u(j))(V • uo)+2S№elk)eiko]dт ( )

T

(j = 2) (

The expressions for the displacement vectors u(j=1) and u(j=2) specified by boundary conditions (15a) and (15b) are given, for example, in [Landau and Lifshits, 1963]. The related components of the stress tensor can be written in the following convenient form:

x

x

1

a(1)

a xx

ao

(1)

a,

yy

ao

(1)

axy

ao a(1) a X2

ao

a(2)

axx ao (2)

3zx2

2nR5

3zy2

2nR5

3zxy

+

+

M

2n(A + m) dy

M d

R(R + z)

2nR5

3xz2

2n(A + m) dx y R(R + z) (ln(z + R))

dxdy

(1)

2nR5 '

“yz

ao

3yz2

(1)

2nR5 ’

a

ao

3z3 2nR5 '

3x

3

a

a

yy

ao

(2)

xy

ao

(2)

2nR5

3xy2

2nR5

3x2y

2nR5

3x2z

+

+

+

M

2n(A + m) M

2n(A + m) M

x d2

R3 + dx2

x d2

R3 + dy2

R + z

R + z

2n(A + m) dxdy

x

(2)

3xyz

a

R + z

(2)

ao

2nR5 ’

ao

2nR5 ’

ao

3x2z

2nR5

dy

R3

2

dy

R5

4

3-r4 ’

y2dy

R5

2

3r2

where

(x - xo)2 + z2

we obtain

a« u XX 2x2z a(2) u XX 2x3

ao nr4 ’ ao nr4 ’

a22) 2z3 a(2) a 22 2z2x

ao nr4 ’ ao nr4 ’

a(1) a xz 2z2x a(2) a xz 2x2z

ao nr4 ’ ao nr4 ’

a(l) ayy Az a(2) ayy Ax

ao n (A + M)r2 , ao n(A + m)'

duZ1) ,, 2

dxo

A ln r2

+B

(x - xo)2

e(X) = C arctan -

x - xo

+B

(x - xo)z

where

A-

1

2^m(A + m)

x V • uo(m^A - —^M - A^M —A^Mexxo

„ 5m I „ A

B =-----------V • uo — + exxo

(51)

C=

nM

1

(52c)

|v • uo(2“^M - ^Aj-^Mexxoj

Thus, the complete solution of the problem under consideration reduces to the calculation of integrals (50a) and (50b).

In the 2-D case (when elastic moduli do not depend on y), simple expressions for tidal tilts and strains can be obtained from formulas (50a, 50b). Integrating components (51) over y from -to to to and taking into account that

These expressions allow one to obtain relatively simple analytical solutions of the problem of tidal deformations in an elastic half-space with heterogeneities independent of y. For example, if inhomogeneities of elastic moduli are specified by a value (5A,5m) in a rectangular parallelepiped infinite along the y axis (x1 < x < x2, -to < y < to, z1 < z < z2), the solution has the form [Molodensky, 1983a]

(52a)

(52b)

n(A + m)

exxo = d(uo)x/dxo ,

and the symbol || indicates the result of double substitution from x1 to x2 in x and from z1 to z2 in z.

6. Thermoelastic Deformation

One of the main shortcomings of the tidal method is its low immunity to the noise produced by local disturbances in the atmospheric pressure, groundwater level, and temperature. The first two effects are either irregular or seasonal and can be eliminated by increasing the lengths of observation series analyzed. On the other hand, the thermoelastic stresses involve strictly periodic diurnal components coinciding in frequency with main solar waves, and their elimination by averaging is ineffective.

We should note that, due to the amplitude modulation of the diurnal temperature variations by the annual wave, this spectrum additionally contains combination frequencies of the diurnal and annual waves with rather high amplitudes. Periods of these waves also exactly coincide with periods of some main lunisolar tidal components (e.g. the Km and Ks waves, whose period is equal to the sidereal day), so that their elimination by methods of spectral analysis alone is also ineffective.

Below we present results of spectral analysis of average temperature variations in the range of tidal frequencies under the natural assumption that, when averaged over a sufficiently long time interval, the diurnal temperature variations are determined by the value of the solar flux incident on the unit area of the Earth’s surface and by the effects of heat capacity of the medium leading to a constant (season-independent) phase delay of temperature. Based on these calculations, we present model corrections to the amplitudes and phases of main tidal waves recorded at the tidal strain-meter Protvino station. Comparison of model results with observations is shown to be effective for obtaining reliable estimates of the thermal expansion coefficient and bulk modulus in the anomalous zone adjacent to the Protvino station.

6.1. Spectral Decomposition of the Solar Thermal Flux.

The value of the solar thermal flux per unit area of the Earth’s surface is determined by the well-known expression

y

x

2

x

x

y2

2

r

2

r

z

2

Table 1. Numerical values of the coefficients in the expansion of the solar thermal flux (6)

L 0 1 -1 2 -2 0 1

N 1 1 1 1 2 2 2

cln 0 0.479 -2.07x10-2 0 0 0 0

L -1 2 -2 0 1 -1 2

N 2 2 2 3 3 3 3

cln 0 9.67x10-2 -5.8x 10-3 0 0 0 0

L -2 0 1 -1 2 -2

cln 0 2.25x 10-2 0 0 2.37x10-2 -5.76x10-3

So cos a at cos a > 0

0 at cos a < 0 , (53)

where So is the heat flow per unit area normal to the direction toward the Sun and a is the angle between the direc-

tion toward the Sun s and the normal to an element of the Earth’s surface. Neglecting the effects of the Earth’s orbit eccentricity, the motion of the vectors n and s in space can be described as

n = i cos(0 - wt) sin $

+ j sin(0 - wt) sin $ + k cos $ , (54)

s = cos fit(i cos £ + k sin e) + j sin fit

and, accordingly,

cos a = (n, s) = sin $ cos(fit) cos e cos(0 - wt)

+ cos(fit) sin e cos $ + sin(fit) sin $ sin(0 - wt) .

Here, i, j, k is the right-handed system of the unit vectors (the unit vector i lies in the equatorial plane, and j is directed toward the vernal equinox point); $ and 0 are the colatitude and longitude, respectively; fi and w are the angular frequencies of the orbital and diurnal rotation of the Earth, respectively; and e = 23.5° is the inclination angle of the equatorial plane with respect to ecliptic.

In the most general case, the expansion of flow (47) in spherical functions can be represented in the form

OO OO n

S = So E E E P"1 (cos $)

1 = o n = o m = -n x (aj(c) cos(m0) cos((mw + 1fi)t)

(s) (55)

+ ainm sin(m0) cos((mw + 1fi)t)

+ b1]C)m cos(m0) sin((mw + 1fi)t)

+ bnl sin(m0) sin((mw + 1fi)t^ .

Because we are further interested only in near-diurnal components of thermoelastic waves, we can set m = 1. The problem in question being fully symmetrical about the Earth’s rotation axis, the coefficients in (49) possess the following properties at m = 1:

a(s) = b(c) ln m ln m ’ (56)

«(c) = b(s) ( )

ln m ln m .

In the case of main diurnal waves, setting m = 1 in (55) and taking (56) into account, (49) can be represented in the following simple form:

OO

S = soEE cln sin(n$) cos(0 - (w + 1fi)t) . (57)

1=o n=o

Numerical values of the coefficients cn in this relation are presented in the Table 1.

Taking into account relation (57) and the fact that the sidereal and solar days Tsid and Tsol are interrelated through the formula

Tsid = Ts(1 - fi/w) , (58)

it is easy to show that, at l = 0, the coefficients cln in (51) describe thermal waves having a period exactly coinciding with the solar day; the l = 1 coefficients describe waves with a period equal to the sidereal day, and the l = 2 coefficients correspond to the frequency which is symmetrical to the frequency of the solar day with respect to the frequency of the sidereal day (i.e. the period of this wave is equal to the difference of two sidereal day and one solar day). All these frequencies are also present in the spectrum of tide-generating forces and can therefore be eliminated from the results of tidal observations only by numerical modelling of thermoelastic strains. Below, we consider several examples of possible effects.

6.2. Equations of Thermoelastic Strains and Their Approximate Solutions

General equations describing the thermoelastic strains in a radially heterogeneous self-gravitating model of the Earth with a hydrostatic distribution of initial stresses can be written as

Li(u, V) = d(kaT)/dxi , (59)

where the operator L is defined by relation (8). Equation

1.04 1.02 1.00 0.98 0.96 0.94

0.92---------------------------------------------------------------------------------------------

0.90---------------------------------------------------------------------------------------------

0.88---------------------------------------------------------------------------------------------

Figure 1. Example of results derived from the finite-difference modelling of the tidal strain effect of a local heterogeneity. (See explanations in the text.)

(59) should be complemented with the boundary condition

(aifcnk - (fcat)ni)s = 0

(60)

(nk are components of the outer normal), describing the absence of stresses at the outer surface of the Earth s.

In the case of local thermoelastic deformations (with the characteristic horizontal size of the heated layer being much smaller than the radius of the Earth), the effects of gravitational forces and initial hydrostatic stresses on the thermoelastic deformations in the Earth is small, and (59) can be replaced by the equation

da ifc dxfc

= d(fcat)/dxi

The effects of local thermoelastic deformations can be significant if the medium contains inclusions with anomalous values of elastic moduli, the thermal expansion coefficient and temperature. Below, we assume that variations in elastic moduli are small compared to their average values in the surrounding medium, whereas no limitations will be imposed on spatial local variations in temperature and the thermal expansion coefficient.

In solving boundary problem (60)-(61) for a nongravitating medium, we may set D = V(J) = D(j) = V = 0. The substitution of (60)-(61) into (15) yields

(61)

which can be considered in the approximation of a homogeneous elastic half-space.

aoux(xo, 0)=/kaT(u(j 1), n)ds

s

- J (u(j=1), V(fcaT))dr = J kaTV • u(j=1)dr

(62a)

Figure 2. Corrections to the x components of tidal tilts and linear strains for an inhomogeneity specified in a rectangular parallelepiped with a 1:5:5 ratio of the sides along the x, y and z axes. (See explanations in the text.)

№“2^ 0)=/kaT <“t'=2,'n)ds

s

- J <u(j=2), V(kaT))dr = J kaTV • u(j=2)dr .

the ratio of the thermal conductivity of the medium to its heat conduction, we have

(62b) z2 = 03w)-1/2 , (66)

where, as before, w = 2n/To and To is the sidereal day.

The divergences of the vectors u(j) in these relations are determined by simple formulas:

ao x - xo

V • u(j=1) =

n(A + m) r2

,(j=2) = ao _z

(63)

u(

r(A + m)

where r2 = (x - xo)2 + z2.

Integrals (62a) are easily calculated analytically for certain simplest models of temperature distribution. For example, if the values kaT are described by a piecewise-continuous function equal to a constant Co = koaoTo in a rectangular volume with vertexes at the points xi,zk (i,k = 1, 2) and is zero outside this rectangle, we have

c (xo, 0) = J kaT V • u(j=1)dr

Co

- (z(lnr - 1)

n(A + m)

+ (x - xo) arctan(z/(x - xo))) |

”■^0) = /kaT V'u"=1’dr

(64)

Co

((x - xo)(lnr - 1)

n(A + m)

- z arctan(z/(x - xo))) || ,

where the symbol | indicates the result of the double substitution from x1 to x2 in x and from z1 to z2 in z.

In calculating the volume integrals (62a, 62b), one should take into account that, in the case of diurnal temperature variations, the depth of the layer of appreciable temperature variations is small compared to the characteristic horizontal size of the integration domain. Expanding (64) in powers of the small ratios z2/(x1 - xo) and z2/(x2 - xo), setting z1 = 0 and differentiating with respect to xo, we obtain

exx = dux (x, 0)/dx|x

^A + — M^ aToz2 n(A + m)

(65)

x2 - xo x1 - xo

where exx is the horizontal component of strain along the x axis at the point xo.

In order to obtain a simple numerical estimate, the value z2 in (65) can be determined as the thickness of the layer in which the amplitude of the diurnal temperature variations decreases by e times. Using the equation of heat conduction in a homogeneous half-space ,3dT/dt = dT/dz2, where f3 is

7. Numerical Modelling Examples

Figures 1 and 2 plot results of finite-difference modelling of tidal strain and tilts for anomalous elastic moduli within a rectangular parallelepiped (x1 < x < x2,-to < y < to, z1 < z < z2).

Figure 1 exemplifies the finite-difference modelling of the tidal strain effect of a local heterogeneity. The spatial coordinate x is plotted on the horizontal axis, and the ratio of the horizontal component of tidal strain in the heterogeneous medium to that in the homogeneous elastic half-space is plotted on the vertical axis. (The model is characterized by an inhomogeneous distribution of the shear modulus M and a constant Lame coefficient A; the shear modulus inhomogeneity is specified within a rectangular parallelepiped that is infinite along the y axis and has a 2:1 ratio of the sides along the x and z axes; the ratio of distances from the outer surface to the upper and lower faces is 0.4; curve 1 is the strain at the upper face of the parallelepiped and curve 2 is the strain at the half-space surface.)

The method of calculating tidal strains in a heterogeneous medium, described above, is particularly effective in relation to 3-D problems involving a complex configuration of the anomalous zone (in this case, the application of the finite-difference and finite-element approaches provides very approximate and often unstable solutions). Figure 2 plots corrections to the x components of tidal tilts and linear strains for an inhomogeneity specified in a rectangular parallelepiped with a 1:5:5 ratio of the sides along the x, y and z axes; in this model variant, the shear modulus M is constant throughout the half-space, and the A value within the parallelepiped is smaller by 90% (the corrections are normalized in such a way that tilts in a homogeneous half-space correspond to the value 7 = 0.7, i.e. maximum and minimum tilts near the vertical faces of the parallelepiped correspond to the values y = 0.9 and 0.5, respectively).

Figures 3 and 4 plot corrections to the x components of linear strains (Figure 3) and tidal tilts (Figure 4) for an inhomogeneity specified within the unit cube (-0.5 < x < 0.5, -0.5 < y < 0.5, h<z< 1 + h; h = 1/10). The values of the Lame parameters within the cube are the same as in the parallelepiped in Figure 2. Curves 1-5 were calculated on the profiles y = z = 0 (1); y = 0.25, z = 0 (2); y = 0.5, z = 0 (3); y = 0.75, z = 0 (4); and y = 1, z = 0 (5). The only distinction of Figures 5 and 6 from Figures 3 and 4 is that they were constructed with the value h = 1.

As an example of the application of relations (64)-(65), we consider possible values of thermoelastic strains for the Protvino station [Latynina and Boyarsky, 1999]. Analysis of long-term observations at this station indicates the presence of well-pronounced seasonal variations in tidal strain am-

2

r

1

1

0.041------------------------------------------------------------------

0.02-

-0.10-

-0.12--------------------------------------------------------------------------

Figure 3. Corrections to the x component of linear strain for an inhomogeneity specified within the unit cube (-0.5 < x < 0.5, -0.5 < y < 0.5, h<z< 1+ h; h = 1/10). (See explanations in the text.)

plitudes (the K1 wave amplitude in summer is higher than its winter value by about 40%). Such an anomaly can be naturally accounted for by the presence of a thermoelastic strain wave that has a sidereal-day period and an amplitude amounting to about 20% of the K1 amplitude and should coincide in phase with the tidal strain variation in summer. With these relations between amplitudes, frequencies and phases, treatment of observation series a few months long fails to separate frequencies of tidal and thermoelastic waves, but variations in the total amplitude reaching 40% of its average value should be clearly distinguishable.

We estimate local variations in the parameter Co that are consistent with such anomalies. A relatively rough (1:50000) map of geodynamic zones separating blocks of the 6th to 10th orders is only available for the area around the Protvino station; this map resolves details whose linear sizes are no less than a few hundred meters. However, as seen from relations (20), the effects of thermoelastic variations attenuate with distance as 1/L3, implying that the effects of local heterogeneities indistinguishable in the map may play the major role. Substituting the typical value z2 ~ 0.2 m into (20) and taking into account geological evidence on the orientation of

Figure 4. Corrections to the x components of tidal tilts for an inhomogeneity specified within the unit cube (-0.5 < x < 0.5, -0.5 < y < 0.5, h<z< 1+ h; h = 1/10). (See explanations in the text.)

Figure 5. The same as in Figure 3 except for h = 1.

the Protvino station relative to geodynamic zones separating blocks, we adopt the values x1 = 500 m and x2 = 1300 m for the respective distances from the tide-recording station to the near and far boundaries of the anomalous zones of orders 8-9. Then, the thermoelastic strain value is

exx ~ 10 4(aTo) .

With the typical values a ~ (10-5 - 10-4), this strain amounts to 5 x 10-9 - 5 x 10-8, which is comparable, on the order of magnitude, with the observed variations in the amplitude of the tidal wave K1 .

Thus, a considerable summer increase in the K1 wave amplitude at the Protvino station is well consistent with the effects of thermoelastic deformations in the anomalous zone near the station.

Relations (19)-(20) and the expansion coefficients of the average solar thermal flux presented in the table can be used

for numerical modelling of average thermoelastic strains with an accuracy sufficient for applications.

Detailed comparison of results derived from the solution of the study problem by various methods and the interpretation of data on local anomalies of tidal tilts and strains in various regions are the subject of the next paper.

Acknowledgments. This work was supported by the Russian Foundation for Basic Research, project no. 04-05-65117.

References

Landau, L. D., and E. M. Lifshits (1963), Theory of Elasticity (in Russian), Moscow.

Latynina, L. A., and E. A. Boyarsky (1999), Seasonal variations in the lunar tide as a model of an earthquake precursor, Vukanol. Seismol. (in Russian), (in press).

Figure 6. The same as in Figure 4 except for h = 1.

Molodensky, S. M. (1977), The effect of horizontal heterogeneities in the mantle on tidal wave amplitudes, Fiz. Zemli (in Russian),

(2), 3-8.

Molodensky, S. M. (1983a), Local anomalies in the amplitudes and phases of tidal tilts and strains, Fiz. Zemli (in Russian),

(7), 3-9.

Molodensky, S. M. (1983b), Gravity field variations due to elastic deformations of the Earth, Fiz. Zemli (in Russian), (9), 3-21.

Molodensky, S. M. (1984), Tides and nutation of the Earth, 214 pp., Nauka, Moscow.

S. M. Molodensky, Schmidt United Institute of Physics of the Earth, Russian Academy of Sciences, Bol’shaya Gruzinskaya ul. 10, Moscow, 123995 Russia

(Received 22 April 2004)