Lateral inhomogeneities —a computational approach

Tectonophys~cs, 57 (1979) T35-T42 o Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands

T35

Letter Section Lateral ~om~eneities

- a compu~tion~

approach

MARKUS BATH Seismological Institute, Box 517, S-751 20 Uppsala (Siueden) (Received March 20,1979;

accepted May l&1979)

ABSTRACT B&h, M., 1979. Lateral inhomogeneities 57: T35--1’42.

- a computational

approach. Tectonophysics,

A computations procedure is developed for the investigation of lateral velocity variations, in an effort to reconcile the geologist’s request for detailed structural information and the seismologist’s need for computationally tractable structures. Based upon great circle arc geometry and including lateral refraction, the methods are applicable to any guided wave propagation, such as of regional crustal waves or of surface waves over any distance.

INTRODUCTION Investigations of the earth’s structure are basically different in geology and seismology. While the geologist generally describes his results in terms of lateral structural variations and co~esponding pe~ologic~ composition, the seismologist is mainly concerned with vertical variations, primarily of wave velocities. For the seismologist, lateral variations generally enter as a second-order approximation. This holds true even on the regional scale, where the seismologist frequently calculates regional events by means of one and the same structural model for his whole region. Even though this is computationally sufficiently accurate in many regions, it often appears inexcusably simple for the geologist. The purpose of the present paper is to contribute to bridging the gap between the geological and the seismological approach by devising a method to investigate lateral velocity variations, still keeping the possibility to treat the structure computationally. The more details one wants to include, the less tractable the s&ucture becomes from the deterministic (non-statistical) compu~tion~ point of view. The restriction to horizontal variations makes the method more applicable to regional and global dimensions than to the purely local dimensions encountered in seismic prospection. The seismological literature contains numerous reports on both global and regional lateral structural variations (see, e.g., B&h, 1959, 1962, Santa, 19650, 1966, and a review by Bgth, 1965). In addition, several methods

T36

have been developed to deal with such geophysical complexities on the global, regional or local scale (see, e.g., Drake, 1972, and Bolt and Smith, 1976). GREAT CIRCLE ARC GEOMETRY

The equation of a great circle arc has the following general form, pressed in the latitude cp and the longitude h of any point on the arc:

ex-

tancp=C1

(1)

sinA--&

cash

C1 and Cz are constants with certain numerical values for any given arc. For their determination, we consider two fundamental cases: (1) The arc is defined by two points with given coordinates, e.g., A(cp,, h,) and E(q,,, A,) in Fig. la. Then we have the correspondence to the “twopoint formula” in plane analytic geometry: CI =C*(&l,cPe,&a

+9o”,xe

+90 tan (Pi sin Xe - tan pe sin X0

C* =G(&z,cPe,L,L?)= (2) The arc is defined

sin&

(2)

-A,)

by one point C(q, , A,) and the angle (Y, measured

Fig. 1. Different arrangements of great circle arcs, discussed in the text. (a) and (b) show the wave propagation along the great circle arcs AE and AH, respectively, neglecting lateral refraction. (c) and (d) show the corresponding wave propagation when lateral refraction is included.

T37

counterclockwise from the meridian (Fig. la). This case corresponds “one-point formula” in plane analytic geometry: C1 =C*((Pc,hc

to the

+90”,(y)

C* =C2(Vc,&,a)=

I

cos A, cos (PC

tan h, tan ff

1

-

(3)

sin (pc I

The eqs. l-3 are generally not given in textbooks on spherical trigonometry or analytic geometry. They can be derived by a straightforward application of the sine and cosine theorems of spherical trigonometry. As special cases, they yield the equation of the equator (V = O”, X = 0” - ?180”) or of any meridian (cp = 0” - *90”, h = constant). By means of the equations given we are in a position to solve any geometrical problem that may enter into Fig. la, of which the following two will be useful in our applications. Problem 1. Calculate the coordinates of the intersection point between two circular arcs. Given: Two circular arcs AE and BD in Fig. la in terms of eq. 1, the constants C1, Cz referring to AE and C’, , Ck referring to BD. Wanted: The coordinates of the intersection point C between AE and BD. Solution: We get A, by eliminating cp between the two equations for AE and BD: tan A, = after

cz -c;

(4)

c, -c;

which (pc is obtained

by the substitution

of X, into the equation

for

AE or BD. Problem 2. Calculate

any distance in Fig. la. Given: Two circular arcs AE and BD in Fig. la in terms of eq. 1 and the coordinates of A. Wanted: The distance AC from A to the intersection point C. Solution: With the coordinates of C obtained as in Problem 1, the distance AC (=A,,) is calculated by the cosine theorem applied to the triangle ACN (IV = north pole): cos Aac = sin pa sin cpc + cos p. cos qc cos(h, In order to make our mathematical duce a function P defined as follows: p(AB)

=

‘OSAa’COSAbe sin Aac

- X,)

expressions

(5) more compact,

we intro-

- cos Aab

(6)

In applications to the earth, the latitudes cp are consistently geocentric, obtainable from the corresponding geographic latitudes (p by the relation: tan cp = 0.993277

tan (p

(7)

T38 WAVE PROPAGATION

Great circle arc propagation In the first approach, we assume that the wave propagates along one great circle arc from the source A to the station 11 (Fig. lb). The equation of AI-I is given as well as the equations for the two discontinuities BD and EG, also assumed to be great circle arcs on the earth’s surface. From the coordinates of A, C, F and Ii, we get the distances AC, CF and FH, travelled in the three sections l-3. The wave velocity is assumed to be constant within each section but to vary from section to section. The travel time equation then reads as follows (t = travel time, u = velocity, a = time term; the subscripts 1, 2, 3 refer to the three sections in Fig. lb): t=

L!E +def Ul

4.

+afh

+a,+a

3

(8)

v3

Equation 8 can be considered as a generalization of the time term method (Bgth, 1978) by allowing for lateral velocity variation in the marker bed. Applied to a number of different paths involving the sections l-3, eq. 8 is solved for the velocities ul, u2 and u3 (as well as for the time terms a, and a3 ) by a least squares procedure. This solution provides a first approximation to our problem, i.e., the three velocities that exhibit the lateral variation. Lateral refraction: one discontinuity If the velocities are significantly different in the different sections, it will be necessary to include the lateral refraction. We consider first two sections (Fig. lc). Using the velocities u1 and u2 found for the straight line propagation (Fig. la) and applying the refraction law and the cosine theorem, we find the following condition for the ray ACE in Fig. lc (il and i ‘1 = incidence angles at the discontinuity between sections 1 and 2): -Ul u2

sin il =-=--= sin i :

P(AB) WW

F(& 1

(9)

This expression is a function only of the longitude h, of the point C. The most convenient way to solve eq. 9 is by a search procedure, i.e., a series of Xc-values are tried until the equation is satisfied. Clearly, there is only one unique solution. Once X, is known, (pc is obtained from the known equation for BD. The distances AC and CE are then calculated by the cosine theorem. After this, a second solution of the travel time equation (8) is obtained, using the revised distances. This solution yields somewhat revised velocity values. With these the solution of eq. 9 is repeated. This iteration is continued until the velocities do not change any more or stay within acceptably small limits. In short, the iteration is this: u from eq. 8 + A from eq. 9 + u from eq. 8 + and so on.

T-39

Lateral refraction: two d~sc~ntinu~ties With more than two sections (Fig. Id), the calculation gets more complicated, but still it can be performed by a search procedure. The longitude X, of C is varied until the refraction law is fulfilled at all discontinuities. The calculation proceeds in the following steps: (1) Assume Xc and get the corresponding ipc from the given equation for Bff. law at C and (2) Calculate ii from the following equation (refraction cosine theorem): sini:

=-

Ul

sini

=-

uz

01

WB)

V2

sin Abe

(IO)

(3) Derive the equation for the great circle arc CF, passing through the point C(P,, h,) and forming an angle = 180” + p + ii with the polar meridian CN, where p is obtained from the sine theorem: cos p =

Cos (Pb sin&

With reference

- xb)

(11)

sin Abe to eq. 3, we then have:

tan a! = tan(l80”

+ p + ii f = tan@ + ii )

(12)

(4) Get the coordinates of F as the intersection point of CF (known from 3) and EG (given). (5) Vary Xe in a search procedure until the following expression for the function G(X,) is fulfilled (refraction law at F and cosine theorem): sin i2 -V2 z = sin ii v3

~P(C=Q) = G@,) WE)

(13)

(6) With all coordinates now known, calculate the distances AC, CF and FH. (7) With these revised distances, repeat the solution of the travel time equation 8. With the new velocities then obtained, the calculation steps (l)-(6) are repeated. In short, the iteration can be described by the following sequence of calculations: v from eq. 8 + A from calculation steps (l)(6) + u from eq. 8 * and so on, oscillating between these two calculations, until the velocities remain constant or stay within acceptably small limits. The essence of the refraction method is to follow the ray path from the source, allowing for refraction at each discontinuity, until the last discontinuity is reached. From here, we join the receiver, and vary hc in a search procedure until the refraction law is fulfilled also at the last discontinuity.

T40

By adding the computational procedure can be immediately between the source A and the case with only one discontinuity

steps (2)-(4) for each discontinuity, the extended to any number of discontinuities receiver 1-I. On the other hand, the simple reduces to the steps (1) and (5)-(7).

APPLICATION

The method offers a device to prove the existence of lateral variations, but it can also be used to disprove suspected variations. An instructive example is provided by explosion investigations of the Swedish crust (B&h, 1971). It was suspected that some velocity difference would exist, especially for Pn, between the eastern lowland and the western highland of Sweden, Therefore, introducing the great circle arc passing through the points 68”N, 24”E and 54”N, 10”E as a demarcation line and applying our calculation method, we arrived at the results listed in Table I. The suspected difference for Pn is clearly denied by the result of this calculation. For Pgl, there is a greater difference between the eastern lowland and the western highland; still the difference is hardly significant. The conclusion is that with the data used, practically no significant lateral velocity variation can be proved for the Swedish area. This is also in accord with B&h (1979), where the same travel times are found useful for the whole of Sweden and adjacent areas in calculating regional events. For example, the travel times for Pgl used by Bgth (1979) deviate by less than 1 set from the two solutions (lowland and highland) in Table I over the distance range from about 100 km to 700 km. On the other hand, it would clearly be incorrect to deny the existence of any lateral velocity variations at all within the Swedish area. But the accuracy and abundance of the data used here are insufficient to demonstrate them beyond doubt. Obviously, we are concerned with scales (cf. Bgth, 1965, p. 491). In order to be re-

TABLE I Examination

of some data on Swedish crustal structure (Bath, 1971)

Wave

Region

Velocity km/set

Pn

total average lowland highland

7.88 f 0.05 7.87 7.90

5.92 * 0.46 5.85 6.01

Pgl

total average lowland highland

6.25 f 0.08 6.05 6.32

0.42 + 0.85 -1.85 1.14

Time intercept’ set

1In terms of eq. 8, these time intercepts are = Za, applicable when both source and receiver are located in the same region. Otherwise, half the values should be used.

T41

vealed, minor lateral variations require denser and more accurate networks (preferably geophone networks) for recording of closely spaced, controlled explosions. If we consider the importance of the scales involved, we are in a better position to reconcile the different approaches of geologists and seismologists, mentioned in the Introduction. Whereas the geologist aims at finding as many structural details as possible, the seismologist is content with a structure that will meet the accuracy of his station network. DISCLJSSION We

summarize our discussion in a number of points, as follows. (I) The developed method can be applied to any guided wave, such as regional crustal waves, teleseismic Pn and Sn, all kinds of surface waves and channel waves. In the case of dispersive waves, the group velocity should be used in calculating travel times, eq. 8, while the phase velocity should be used in calculating refraction, eq. 9 and following (Bullen, 1963, p. 107, and Bikh, 1968, p. 48). Therefore, for dispersive waves, the method should be supplemented by relations between group and phase velocities. However, as velocities enter only in the form of ratios in the refraction formulas, it is possible to use the group velocities also here with acceptable accuracy. (2) A salient feature of the method is that structures are derived in a form immediately suitable for further computation. Hence, a derived, laterally inhomogeneous structure can be applied for a more accurate location of regional events by an iteration procedure (cf. Herrin and Taggart, 1962). In the case of regional events, the same waves are used in the location calculation as in our method, and therefore a more efficient way to derive a structure would be by the use of controlled explosions. For the surface waves from teleseisms, this complication does not arise, as the location is then based on the independent P waves, and natural events will then be sufficiently accurate. (3) The method may be applied to investigations of the lateral variation of other distance-dependent properties, e.g., the attenuation of guided seismic waves. (4) The calculation procedure offers an advantage over direct measurements, e.g., on maps or globes, by being more accurate. All calculations involved can be conveniently handled by any computer of a sufficient capacity. Moreover, the use of spherical trigonometry makes the method applicable to any scale, both regional and global. (5) Geometrically, the restriction of discontinuities (structural boundaries) to great circle arcs could be considered as a certain disadvantage. However, this is circumvented by composing the boundary of a sequence of great circle arcs, each of which can be made short enough to suit any given boundary. Any number of boundaries of any shape will make a mosaic-pattern approach possible. Even though some boundaries could be more easily represented by other contours, e.g., small circle arcs, their use would not change

T42

the methodological principles. For very complicated boundaries, it may still be more convenient to combine our method with readings from maps. (6) Geophysically, the location and the shape of the boundaries are not given a priori, but they have to be chosen from geological, topographical and other features. Alternatively, a regular net of boundaries may be a useful, first approach. (7) We have assumed sharp discontinuities. However, in nature, lateral discontinuities may more often be gradual than abrupt. Our assumption is no doubt a simplification, but it will still demonstrate the main features. REFERENCES BBth, M., 1959. Seismic surface wave dispersion: A world-wide survey. Geofis. Pura Appl., 43: 131-147. 1962. Crustai structure in Iceland and surrounding ocean. ICSU Rev., B&h, M., 4: 127-133. Lateral inhomogeneities of the upper mantle. Tectonophysics, B&h, M., 1965. 2: 483-514. Bgth, M., 1968. Mathematical Aspects of Seismology. Elsevier, Amsterdam, 415 pp. B&h, M,, 1971. Average crustal structure of Sweden. Pure Appl. Geophys., 88: 75-91. B&h, M., 1978. An analysis of the time term method in refraction seismology. Tectonophysics, 51: 155-169. B&h, M., 1979. Earthquakes in Sweden 1951-1976. Swedish Geol. Survey, C 750, Yearbook 72, No.‘12, 79 pp. Bolt, B.A. and Smith, W.D., 1976. Finite-element computation of seismic anomalies for bodies of arbitrary shape. Geophysics, 41: 145-150. Bullen, K.E., 1963. An Introduction to the Theory of Seismology. Cambridge Univ. Press, Cambridge, 381 pp. Drake, L.A., 1972. Love and Rayleigh waves in non-horizontally layered media. Bull. Seismol. Sot. Am., 62: 1241-1258. Herrin, E. and Taggart, J., 1962. Regional variations in Pn velocity and their effect on the location of epicenters. Bull. Seismol. Sot. Am., 52: 1037-1046. Sant6, T., 1965a. Lateral variation of Rayleigh wave dispersion character. Part I. Observational data. Pure Appl. Geophys., 62: 49-66. Santa, T., 1965b. Lateral variation of Rayleigh wave dispersion character. Part II. Eurasia. Pure Appl. Geophys., 62: 67-80. Santa, T., 1966. Lateral variation of Rayleigh wave dispersion character. Part III. Atlantic Ocean, Africa and Indian Ocean. Pure Appl. Geophys., 63: 40-59.