Dissipation of Energy Radiated by Earthquakes

SUMMARY. - Representing radiated seismic energy as a sum of normal modes we are able to develop approximate equations describing the dissipation of this energy in the earth by solid friction. In particular we calculate total energy radiated and spherically averaged dissipation as a function of depth. In a sample (long period) model study we find that energy density as a function of radius and frequency have simple forms and smooth variations with source depth and mechanism as long as the event's spectrum is dominated by fundamental mode energy. Further, we find that shallow focus events are an order of magnitude more energetic than deeper events for the same moment and that dip slip events are up to a factor of six more energetic than strike slip events for the same moment and depth.

SUMMARY. -Representing radiated seismic energy as a sum of normal modes we are able to develop approximate equations describing the dissipation of this energy in the earth by solid friction. In particular we calculate total energy radiated and spherically averaged dissipation as a function of depth. In a sample (long period) model study we find that energy density as a function of radius and frequency have simple forms and smooth variations with source depth and mechanism as long as the event's spectrum is dominated by fundamental mode energy. Further, we find that shallow focus events are an order of magnitude more energetic than deeper events for the same moment and that dip slip events are up to a factor of six more energetic than strike slip events for the same moment and depth.
(*) Institute of Geophysics and Planetary Physics Scripps Institution of Oceanography University of California, San Diego La folia, California 92093. INTRODUCTION. In a recent paper Kanamori (1977) has developed a new magnitude scale, M w , based on scalar moment Mo, minimum strain energy drop Wo, and the Gutenberg-Richter magnitude energy relation.
With the magnitude scale the 1960 Chilean earthquake is assigned M w = 9.5. Since the largest earthquakes are dominant contributors to radiated energy (Stacey, 1977, Ch. 5) it is intriguing that past estimates of the energy release in earthquakes may be an order of magnitude too low.
In this paper we present a method for estimating not only the energy radiated (McCowan and Dziewonski, 1977) but also the radial distribution of the dissipation of the energy radiated. Our method is much more elaborate, and expensive, than the usual methods that involve a logarithmic relationship between energy and moment or magnitude Furthermore, cur method requires that one have a good model of the Earth, including both elastic and anelastic parameter distributions.
GROUNDWORK. Assume a spherically symmetric, self gravitating, non-rotating earth model with density p (/•), and Lamé parameters [i (r) and A (/•). In the notation of Gilbert and Dziewonski (1975) (hereafter references as M), the eigendisplacement of the }"' normal mode is s, (r) such that where w,-2 is the j'" eigenfrequency, S,k is the Kroneker delta and * means complex conjugation. The displacement, u (r, t), due to an external Wo~Mo/(2x 10 4 ) log Wo= 1.5 M,"+ 11.8 where a, = w ; /2Qy is the damping parameter and Q, the quality factor of the /"' mode, H (/) is the Heaviside function, (*) means convolution, and is the source excitation of the /"' mode. where 'Z', (OJ) is the Fourier transform of </*, (/) and A, = OJJ + /«,. We can think of three kinds of singularities of U (r, w) that contribute to u (r, t). The singularities of IP,-(co) represent the deviation of the time behavior of / (r, t) from that of a step function. For our purposes the phase of ' P,-(OJ) is not important and we shall neglect the singularities of x Pj (cu). The pole at w = 0 gives rise to the statical part of u{r,t), the part remaining as f-» <». This part is also unimportant for cur purpose because it does not represent radiated energy. We shall neglect the pole at a> = 0. Finally, we have the poles at w = <r ; -and co=-Oj*. These represent the free oscillations and do account for the radiated energy. The part of a (r, t) arising from these poles is denoted I/» (r, 0 u 0 (r, t)= Re °ll (r, t) [5] °H ('•, t)= I Sj (r) ^ e ia j l H (/).
In other words, there are two approximations involved in transforming [2] to [5], First, the constant part of the mode shape (the static part of (he displacement) has been dropped. Second, we have approximated fjjj (/) as "P, (<7,) S (I) (or x I'j (a>) as a constant, !P; (IT,-)). This is partially justified because ¡'/ / /(w)j is very smooth compared to a normal mode resonance peak (i. e. 0/ (/) must be concentrated at the origin time of the event). But the phase of (a>) may vary rapidly due to source complexities. However, source phase does not enter into the energy calculation. Finally, note that in [5] we have changed our notation to a complex exponential representation. TOTAL DISSIPATED ENERGY.

ENERGY DISSIPATION AS A FUNCTION OF RADIUS.
To determine where in the earth the radiative dissipation takes place we introduce the dissipative stress power density P D (Malvern, 1969, p. 267) which is obtained by contracting the dissipative part of the stress tensor with the strain rate tensor.
In terms of [5] the strain tensor is the real part of E (r, t) = 2 (y 0,1 + A,) % (a,) e'">' [10] where dj=V-SJ(r) and A, is the deviatoric strain for the j' h mode, 1 is the unit tensor. The stress tensor is the real part of T (r, t)=Z (x0, 1 + 2/iAj) Vj (*,) e'"i' [11] i 2 and x = X + -fi is the bulk modulus. We regard x and fi as complex li-' (r)=A£o" 1 (r) (1 -iQ," 1 (r)). [12] It is incorrect to suppose, as we have here, that the model attenuations, Qu and Q x , are frequency independent as this leads to non-causal results. However, the best current evidence (Akopyan, Zharkov) and Lyubimov (1975);Liu, Anderson and Kanamori, (1976) indicates that Q" and Q x vary slowly (if at all) over the whole seismic band and may be considered constant over the normal mode band (.3 to 16.7 mHz). We feel it consistent with other approximations to assume a frequency independent Q model.

The more common form of [12] is
where relative errors 0 (Q^1) are neglected. Combining [11] and [13] F. GILBERT -R. BULAND the dissipative part of the stress tensor is the real part of 7'° (r, 0=2 (i (x" Qx" 1 0; 1 + 2/Uo <V' 4/) (o>) ' [14] To this degree of approximation the strain rate tensor is the real part of do not require a spherically symmetric earth model. They would be equally valid if the model was a function of r instead of r. However, as a practical (computational) matter there is currntly no advantage in preserving such generality. In the following spherical symmetry is necessary.
We introduce the spherically averaged, time averaged dissipative stress power density We invoke the orthogonality of the spherical harmonics in the representation of Sj (r) to evaluate [17] P°"(r, 0=1/2 T° (r, 0: £* (r,t). [16] [17] where [19] Expressions for A.' , and M, have been given by Backus and Gilbert (1967, Appendix B) in terms of the scalars used to represent Sj (r).
In [18] we have neglected the interaction between overtones with the same angular order. With overtone numbers n and m, and using the representation [2] temporarily, the interaction involves terms such as The average over a few cycles of the cosine term with argument (w,! + w,") t is zero and can be neglected. The same is true of the term with argument (ou" -w,") t unless \a)" -cj m \ <a" + a m . In that case we use the fact that Two overtone frequencies can be arbitrarily close if the modes are « disjoint », i. e. each mode has an eigenfunction that is significantly different from zero only where the eigenfunction of the other mode is nearly zero. Thus, the contribution of the interaction terms to r 3 )D in [18] can justifiably be neglected, a fact that is not self-evident when the complex exponential representation is used.

[22]
Note that O, is the quality factor of the /'"' normal mode while Q,, (r) and Q x (r) represent the attenuative properties associated with the shear and bulk moduli respectively at radius r.
The utility of the expression [22] is that one can calculate not only the total dissipated energy ¿° but also its radial distribution r £°. That is, one can calculate how much energy is dissipated in different regions such as the inner core, lower mantle, crust, etc. Furthermore, one can quickly investigate the effect of the epicentral depth and mechanism of the source on the radial distribution of dissipated energy. As an example of the method, we have performed a limited model study for the effects of varying depth and mechanism of the source on seismic energy dissipation. Although parameter space was not searched exhaustively the results give a good qualitative feel for what is happening. Indeed the results are not only more informative but also are more general than one might suppose for the following reasons. First, of all physically realizable seismic source mechanisms only a handful are at all common. Second, variations in the dip, slip, and depth of a source are found to cause smooth changes in the energy release. Third, since the source phase does not enter into the calculation there is a good deal of symmetry in the source parameters. Identical pictures are generated, for instance, by any strike, by interchanging the slip and nodal planes, or even reversing the slip direction. And finally, any geographical location will give the same results since [22] is spherically averaged.
We began with earth model 1066A (M, Table 5) and a simple Q model (M, p. 199). We calculated all elastic-gravitational free oscillations with periods longer than 66 s. Eigenfunctions and frequencies of a purely elastic model were used as the effects of physical dispersion are on the order of other approximations made above. The sources were all step functions in time and delta functions in space (point sources). Each had a moment of \0 20 N'm. Normal mode excitations were calculated as in M (Section 2) for interior sources and as in Appendix A for shallow (surface) sources.
The results of our experiments are presented in Fig. 1 -8. Each figure consists of two pictures, energy density as a function of radius (in units of a, the radius of the earth) and energy density as a function of frequency (in millihertz) running mean averaged over 1.67 mHz windows. Total energy, at periods longer than 66 s., has been normalized to unity. The source used to generate each figure is designated by a shorthand four-vector: (depth in kilometers, dip in degrees, slip in degrees, and total radiated long period energy in Joules). Dip has the usual definition. Slip describes the direction of the slip on the fault plane. It is defined so that 0° is pure strike slip, +90° is normal faulting, and -90° is thrust faulting.
Examination of the eight plots of energy density versus radius shows that the zeroth order shape of the curves is controlled by the Q model. Indeed, these curves are linear in (r). Thus the low model Q" (=115) above .9 results in high attenuation in the upper mantle. The intermediate Q"( = 400) above r~ .55 results in small losses in the lower mantle. And the high Q^(=1000) in the inner core translates to negligible losses. Superimposed on Q^" 1 is a first order structure that is a function of source depth alone. Thus all plots of shallow sources (Figs. 1-5) are very similar. They have a sharp peak near the surface, a broader peak at r~ .97, and a roughly exponential fall off to the core (of course this tail is amplified about 4: 1 in the upper 650 km. by the Q" contrast). The high amplitude spike is due to shallow fundamental mode surface wave energy amplified in the low velocity crust. The broad peak corresponds to surface wave energy trapped in the low velocity zone (LVZ). The fall off is due to the exponential tails of the well excited upper mantle surface wave modes. That is, fundamental mode energy dominates completely for shallow sources.
At intermediate depth the picture is the same (Figs. 6 and 7) with two exceptions. First, the surface spike has been rduced relative to the LVZ peak. Second more energy has been dissipated below the LVZ. For deep sources the picture (Fig. 8) is different indeed. The surface spike is gone as the deep source is below the turning points of the modes that composed it. The LVZ peak is still visible, but is now comparable to higher mode energy deeper in the mantle. The deeper the source the more the higher modes are excited. For very deep sources there is as much energy in modes as in the very long period fundamental modes.
The companion plots of energy density versus frequency show considerably more variation with source mechanism. However, all shallow and intermediate sources are characterized by curves that vary smoothly and apparetnly predictably with frequency except at the longest periods. This is a result of the dominance of the fundamental modes. For a shallow source one would expect the fundamental mode excitations to vary smoothly along the branch to an arbitrarily high frequency (i. e. the given curves may be extrapolated). For intermediate depth sources the character of the curv must change dramatically at some higher frequency where the equivalent ray turning point becomes structure due to the comparable excitation of higher modes and long shallower than the source depth. The deep source shows a more complex period fundamental modes.
Figs. 1-3 demonstrate the variation in energy density versus frequency with dip for purely dip slip sources. Between Figs. 1 and 2 (Fig. 1 represents a dip of 67.5° as well as 22.5° as both have the same absolute excitation) we see that dip effects the results very little until one approaches vertical dip slip (Fig. 3). Vertical strike slip (Fig. 4) has a very different structure because of a different frequency dependence of the fundamental mode excitations. When vertical strike slip and vertical dip slip are mixed (Fig. 5) the result, as one would expect, is a combination of Figs. 3 and 4. The intermediate depth pictures (Figs. 6 and 7) support the finding that dip in a dip slip event is relatively uniportant.
The final parameter we consider in our model study is total seismic energy released (at periods longer than 66 s.). Again the first order variation is a depth dependence. Intermediate depth sources are almost five times as energetic as deep sources of the same moment. Shallow sources are an order of magnitude more energetic than intermediate dpth ones. Further, for shallow sources there is a large variation among various mechanisms. Thus vertical dip slip is only half as energetic as a 45° dip slip source and is three times as energetic as a strike slip source.
As a further example, we have applied the calculations described above to the 31 July 1970 Colombian earthquake. We used normal modes with periods longer than 80 s. and the source function of M ( Fig. 6). The result is displayed in Fig. 9. Notice the similarity between Fig. 8 and Fig. 9. The additional complexity in the source has changed the result very little except to multiply the energy density versus frequency curve by different absolute squared source spectrum.

DISCUSSION.
We have presented a method of calculating where in the earth radiated seismic energy is dissipated by solid friction. Our equations are correct tc 0 (Q~'). The model calculations, presented above, reveal a pattern which may be used to advantage in the more speculative calculation of a seismic component to the upper mantle heat budget.
It is well known (Stacey, 1977, Ch. 5) that the annual seismic energy release is dominated by a few of the very largest events. We have shown additionally that one may neglect all but shallow events and all but the greatest shallow strike slip events. Further, we have shown that the shape of the energy dissipation as a function of radius (at least at long periods) is relatively independent of source mechanism.
By making a few assumptions it is quite possible that one could reasonably estimate the average annual seismic energy dissipated at each depth in the mantle. One would need to estimate source mechanism and absolute source spectrum for a handful of events. Notice that seismic efficiency does not enter into the calculation if the moment density tensor is estimated from radiated seismic energy. If one could estimate the contribution of wavelenghs on the order of the source dimension and shorter the job would be done.

ACKNOWLEDGEMENT.
This research has been sponsored by the National Science Foundation under grant NFS EAR77-00897.
APPENDIX A Seismic sources where the sources dimension is on the order of the source depth (that is almost any large shallow earthquake) cannot reasonably be represented as points even for wavelengths long compared with the source dimension. To see this we consider the modal excitation of a point source. Following M (Section 2.1) we rewrite [3] as where M is a moment density rate tensor such that and A is a Dirac delta function 4 <r)=* (r, r")= ^^ [28] for a source at r 0 = (r 0 , 0 O , (po is the strain due to the /"' mode.
For an interior source the volume integral in [29] contributes at the source and the surface integral vanishes. But at the free surface certain strain components vanish in order to satisfy the traction free boundary condition. Hence, for a shallow point source corresponding components of M excite neglibible seismic energy. For example, a near surface vertical dip slip equivalent point source would excite no modes. Allowing the source to be at the surface is no solution as the surface integral in [29] is then infinite.
Clearly, we need to consider a distributed source. A proper three dimensional equivalent body force density distribution is beyond the scope of approximations made above. Rather we relax our definition of A (r) in the radial dimension only. oL [31] r sin tj where ' F (r) dr -1 [32] so that seismic moment is preserved and a J'F(r)dr=0 [33] u so that momentum is conserved. F (r) is considered to be delta-like so that the source is concentrated about the mean source depth. Now for a shallow source, r 0~( a,do, 0o), [29]  -1/2 F (a) S/* (r 0 ) + (r 0 ) r): M (r 0 , 0 From [52] and the fact that F (r) is delta-like we see that F" 1 (r) has the dimensions of length. Without considering the actual shape of F (r) we have taken F" 1 (a) in [34] to be 50 km. for the « surface » source used in the calculations.  GILBERT, F., and DZIEWONSKI, A. M., 1975. -An application of normal mode theory to the retrieval of structural parameters and source mechanism from seismic spectra. Phil. Trans. Roy. Soc.,