S. M. Molodensky
Schmidt United Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia
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 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 d (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.
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 p_{0} is compensated for by gravity:
(1) |
where r is density; W is the unperturbed gravitational potential, satisfying the Poisson equation
(2) |
and G is the gravitational constant.
Elastic tidal deformations displace an element of the medium from a point
determined by the radius vector
(3) |
(4) |
(5) |
where d_{ik} is the Kronecker delta and s_{ik} 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
(7a) |
(7b) |
Applying the curl operator to the left- and right-hand sides of (1) yields
Adding the identically zero term
to the left-hand part of Equation (7a), this equation takes the form
(8) |
(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).
To describe tidal deformations in a radially and laterally heterogeneous medium, we use Equations (8) and (7b) with the inhomogeneous boundary conditions
(9a) |
(9b) |
where
S and
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
(10a) |
(10b) |
where
F (
(10c) |
(10d) |
which satisfy the boundary conditions
(11a) |
(11b) |
(11c) |
(11d) |
Here
(12) |
is the density of the single layer on the Earth's surface, e_{r} is the unit radius vector and y^{m}_{n} are normalized spherical functions.
Taking into account the explicit form of the operator
(13a) |
reduces to the surface integral
(13b) |
where the density of the single layer D_{1} is determined by a relation similar to (12):
(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
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 solution ( u^{(j)}, V^{(j)}) ); the additional terms
determine the difference between the gravitational energies of the single layer D_{1} in the field with the potential R^{(j)} and the single layer D^{(j)} in the field with the potential R_{1}.
If the external body forces
F(
(15) |
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(
Indeed, if external surface forces are absent, we have
s_{1} n = 0. Let
the solutions
( u^{(j)}, V^{(j)} ) be specified
by the boundary conditions
s^{(j)} n =
e^{(j)}d ( r - r_{0})
and
D^{(j)} = 0, where
e^{(j)} is an arbitrarily oriented vector,
d( r - r_{0}) is the 2-D
delta function, and
(16) |
which completely defines the three components of the sought-for vector u_{1} at any point r_{0} of the surface s.
Analogically, specifying the solutions ( u^{(j)}, V^{(j)} ) by the boundary conditions s^{(j)} n = 0 and D^{(j)} = d ( r - r_{0}) and integrating the right-hand part of (15) over angular variables, we obtain the relation
(17) |
which completely defines the potential variations on the Earth's surface.
Tidal displacements and the potential in spherically symmetric models of the Earth are determined by the Love formula
(18) |
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 r_{0}(r) and elastic moduli l_{0}(r) and m_{0}(r) ) is replaced by a spherically nonsymmetrical model (with distributions r_{0}(r) + dr, l_{0}(r) + dl and m_{0}(r) + dm ), the solutions ( u_{0}, R_{0} ) of Equations (8) and (7b) acquire increments ( u_{1}, R_{1} ).
In the first approximation of the perturbation theory, the equations of tidal strain in a laterally and radially heterogeneous medium can be written as
(19a) |
(19b) |
where
(19c) |
and
(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:
(20a) |
(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:
(21a) |
and
(21b) |
Here f_{i}(r) is a system of six functions defined by the ordinary differential equations
(22) |
with the coefficients
(23) |
The symmetry properties of the coefficients a_{ik} 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:
(24) |
Then, substituting (20) and (21) into (13), we find
(25) |
Here
(26) |
These relations completely define all of the sought-for coefficients in expansions (20a) and (20b).
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 x_{nm} = r^{n} y^{m}_{n} (, j). In Cartesian coordinates x_{i}, they satisfy the equation
(27) |
and the Laplace equation
(28) |
Expressing the displacement vectors u_{0} and u^{(j)} through x_{nm}, we obtain
(29a) |
(29b) |
where x_{0} = r^{2} y^{j}_{2} (, j); j = 0, 1 and 2 for long-period, diurnal and semidiurnal waves, respectively; and the functions h^{(j)}_{n}, t^{(j)}_{n} and h^{(0)}, t^{(0)} depend on the radius r alone. Comparing (18), (21a) and (29), these functions can easily be expressed through f^{(j)}_{i}(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
(29c) |
where
(29d) |
Using the Gauss theorem and taking into account the condition of vanishing stresses at the outer surface,
we obtain
(30) |
The substitution of (19d) into (20) yields
(31a) |
where
(31b) |
and
(31c) |
Substituting (29a) and (29b) into (31c), we obtain
(32) |
Taking into account (27), we have
(33) |
The expansion of the lateral inhomogeneities of elastic moduli in series of harmonic functions yields
(33a) |
Taking into account (32)and (33), the integrals in (31b) and (31c) can be represented as the sum of terms containing integrals of products of three harmonic polynomials, scalar products of their gradients, and the convolution of the tensors consisting their second derivatives and having the form
(33b) |
where
(33c) |
For calculating the integrals of products of three harmonic polynomials, we substitute the relation x_{nm} = r^{n} y^{m}_{n} (, j) into the first term in (33b) and obtain
(34) |
Here A_{lpnmn0m0} is the integral of the product of three spherical functions with the respective indexes ( l, p ), ( n, m ) and ( n_{0}, m_{0} ) taken over the surface of a unit sphere.
The second term in formula (33b) contain the product of the homogeneous harmonic polynomial x_{lp} and the scalar product of the polynomial gradients ( x_{nm}, x_{0} ) and can be calculated as follows. Denote
(35) |
Integrating (35) by parts and taking (27) into account, we find
(36) |
and two similar relations obtained from (36) by the cyclic permutation of the polynomials x_{lp}, x_{nm}, x_{n0m0}. Solving these equalities, we find
(37) |
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
can be calculated in a similar way. Denoting
(38) |
and integrating (38) by parts, we obtain
(39) |
where
(40) |
Taking (27) into account, we have
(41) |
and, therefore,
(42) |
Integrating J^{n0m0 (3)}_{lpnm} by parts, we obtain
(43) |
Applying the cyclic permutation of indexes ( l, p)(n, m) (n_{0}, m_{0} ) to this relation and determining J^{n0m0}_{lpnm} from the resulting system of three equations, we find
(44) |
or, substituting (37),
(45) |
The substitution of (42) and (45) into (39) yields
(46) |
With the help of (37) and (36), expression (32) is transformed into the final form
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
(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 ( x_{0}, y_{0}, 0 ):
(49a) |
(49b) |
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
(50a) |
(50b) |
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:
(51) |
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 and taking into account that
where
we obtain
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 ( dl, dm ) in a rectangular parallelepiped infinite along
the
y axis ( x_{1}
(52a) |
(52b) |
where
(52c) |
and the symbol | indicates the result of double substitution from x_{1} to x_{2} in x and from z_{1} to z_{2} in z.
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 K_{m} and K_{s} 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 strainmeter 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.
The value of the solar thermal flux per unit area of the Earth's surface is determined by the well-known expression
(53) |
where
S_{0} is the heat flow per unit area normal to the direction toward
the Sun
and
a is the angle between the direction toward the
Sun
(54) |
and, accordingly,
Here,
i, j, k is the right-handed system of the unit vectors (the unit vector
In the most general case, the expansion of flow (47) in spherical functions can be represented in the form
(55) |
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:
(56) |
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:
(57) |
Numerical values of the coefficients c_{n} in this relation are presented in the Table 1.
Taking into account relation (57) and the fact that the sidereal and solar days T_{ sid} and T_{ sol} are interrelated through the formula
(58) |
it is easy to show that, at l = 0, the coefficients c_{ ln} 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.
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
(59) |
where the operator L is defined by relation (8). Equation (59) should be complemented with the boundary condition
(60) |
( n_{k} 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
(61) |
which can be considered in the approximation of a homogeneous elastic half-space.
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
(62a) |
(62b) |
The divergences of the vectors u^{(j)} in these relations are determined by simple formulas:
(63) |
where r^{2} = (x - x_{0})^{2} + z^{2}.
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 C_{0} = k_{0}a_{0}T_{0} in a rectangular volume with vertexes at the points x_{i}, z_{k} (i, k = 1, 2) and is zero outside this rectangle, we have
(64) |
where the symbol | indicates the result of the double substitution from x_{1} to x_{2} in x and from z_{1} to z_{2} 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 z_{2}/(x_{1} - x_{0}) and z_{2}/(x_{2} - x_{0}), setting z_{1} = 0 and differentiating with respect to x_{0}, we obtain
(65) |
where e_{xx} is the horizontal component of strain along the x axis at the point x_{0}.
In order to obtain a simple numerical estimate, the value z_{2} 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 b T/ t = T/ z^{2}, where b is the ratio of the thermal conductivity of the medium to its heat conduction, we have
(66) |
where, as before, w = 2p/T_{0} and T_{0} is the sidereal day.
Figure 1 |
Figure 2 |
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 l; 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 l 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 g = 0.7, i.e. maximum and minimum tilts near the vertical faces of the parallelepiped correspond to the values g = 0.9 and 0.5, respectively).
Figure 3 |
Figure 4 |
Figure 5 |
Figure 6 |
We estimate local variations in the parameter C_{0} 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/L^{3}, implying that the effects of local heterogeneities indistinguishable in the map may play the major role. Substituting the typical value z_{2}0.2 m into (20) and taking into account geological evidence on the orientation of the Protvino station relative to geodynamic zones separating blocks, we adopt the values x_{1} = 500 m and x_{2} = 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
With the typical values a(10^{-5}-10^{-4} ), this strain amounts to 510^{-9}-510^{-8}, which is comparable, on the order of magnitude, with the observed variations in the amplitude of the tidal wave K_{1}.
Thus, a considerable summer increase in the K_{1} 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.
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).
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.