Geophys. J. R. astr. SOC.(1971) 25, 49-69.

A Linear Model for the Yield Dependence of Magnitudes Measured by a Seismographic Network Ulf Ericsson (Received 1971 February 4)

Summary The connection between the yields of six underground nuclear explosions in the Pahute Mesa in Nevada and the body and surface wave magnitudes measured at 19 Canadian seismographic stations is described by a simple statistical model. It describes the expected magnitudes as linear functions of the logarithm of the yield and attributes the apparently normally distributed magnitude variations between stations to wave path dependent level differences between the transfer functions from the source area to the stations. The model assigns the same scaling exponent for yield to all stations. The exponent is, however, larger for surface wave magnitudes than for body wave magnitudes. The station residuals are partly attributed to differences, between explosions, in the energy coupling in the source area. The residuals are roughly normally distributed and the stations can be divided into groups with high and low residual variances. Applications of the model to the prediction of station magnitudes from given yields, to the calculation of maximum likelihood estimates of yields from measured station magnitudes and to the identification of underground nuclear explosions by simultaneous measurements of body and surface wave magnitudes are briefly discussed.

1. Introduction The connections m ( W ) and M ( W ) between the yields or energies W of underground nuclear explosions and the seismological magnitudes m and M of the generated body and surface waves have several practical applications. One is the prediction of magnitudes from contemplated yields, a matter of political importance in an eventual ban on underground nuclear explosions above some threshold magnitude. Another technical application of some interest is the estimation of yields from observed magnitudes. A third application is to the analysis of simultaneous measurements of m and My a matter important to the possibilities to identify nuclear explosions and earthquakes under an eventual test ban on underground nuclear explosions. In addition the understanding of m ( W ) and M ( W ) should also be helpful in the study of earthquake energies. On experience with M ( W ) nothing seems to have been published so far. There are, however, several papers on experience with m ( W ) , see e.g. Romney (1959), Latter & Martinelli (1959), Pasechnik et al. (1960), Riznichenko (1960), Carpenter, Savill & Wright (1962), Werth, Herbst & Springer (1962), Romney (1962, 1963), Werth & Herbst (1963), Hankins (1964), Green, Frosch & Romney (1965), Springer (1966), 49

D

50

Ulf Ericsson

Hays (1969), Evernden (1969, 1970) and Basham (1970a). From these it appears that m can be expected to be linear or nearly linear in log,, W over at least one unit of log,, W and that there are large numerical variations in and between the results obtained by different authors. Part of these can be attributed to different definitions of m and different instruments, but from recent measurements by fairly homogeneous station networks, such as the WWSSN and LRSM, it has become clear that a significant part of the magnitude variations is due to the inhomogeneities of the Earth's crust and mantle. This can be seen in Table 1 below, where some means and standard deviations of observed explosion magnitudes are collected from various sources. The determination by Marshall et al. for the Long Shot explosion is particularly significant, as the WWSS network seismograms were read by one analyst in a deliberate attempt to be consistent. Contour maps published by Jordan, Black & Bates (1965), by Reagor, Gordon & Jordan (1968) and by Lambert et al. (1969) show the geographical distribution of the m and M or corresponding amplitudes/period for such large explosions. There are large deviations from radial symmetry around the sources and in a fashion resembling meteorologic high and low pressure areas, with trend areas in between. The anomalies appear to be geographically tied but shift with source position. The sometimes strong dependence of regional M-anomalies on the position of the source has also been demonstrated by McGarr (1969), who investigated the distribution of the M obtained from 25 WWSSN stations in the conterminous US and generated by earthquakes in the Pacific and Atlantic. McGarr explained the dependence of the M-anomalies on source position by horizontal refraction at crust and mantle irregularities along the wave path. One particular conclusion from such anomaly pattern changes with changes of the position of the source is that the use of additive source corrections and station corrections is insufficient. A geographical anomaly pattern defined by a set of station corrections is not changed in shape by the use of additive corrections for source position, only in level. It is like changing a depth chart from fathoms to feet to meters, depending on where it is published, and not like making a different chart. From the above it appeared to the author that a realistic or practically useful description or model of m ( W ) and M ( W ) would have to distinguish between various magnitude definitions and station instrumentations, would have to be specific to the position of source and station and would have to consist of a linear or nearly linear deterministic or functional relationship between log W and the expected magnitude plus a stochastic part, describing the magnitude variations because of noise and because of variations between explosions of the coupling of the explosion energy into the ground. The present study is an attempt to obtain such a model. Data useful for such a study would be a collection of magnitudes from a fixed network of specified stations, for several explosions with a wide range of well determined and published yields in each of several epicentre areas, all with published co-ordinates, depths and geologies. The production of such data appears, however, to have been somewhat limited. Well determined, that is official radiochemical, yields have been published only for some 50 specified explosions, all in the US and most in one epicentre area in Nevada. Only about 20 of these had yields above 20 kt and were thus strong enough for widespread teleseismic observation. 2. The yields

The yields given by the US testing agency for the six explosions used here, were taken from a publication by Higgins (1970) and are reproduced in Table 2, together with code names, dates, epicentres and hypocentre depths below surface. The explosions occurred in the Pahute Mesa in the north-western part of the Nevada test

2 3 4 5 6

1

Index

Date 66 05 06 66 06 30 66 12 20 67 05 23 67 05 26 68 04 26

Reference USCGS EDR 89/65 Jensen et al. (1966) Marshall et al. (1966) Lieberman & Pomeroy (1967) USCGS EDR 64/69 Lieberman & Pomeroy (1967) Reagor el al. (1968) Lambert et al. (1969)

Name Chartreuse Halfbeak Greely Scotch Knickerbocker Boxcar

Explosion Long Shot Long Shot Long Shot Saphir Milrow Saphir Gasbuggy Long Shot

N. Latitude 37" 20' 53" 18 57 18 07 16 30 14 53 17 44

Table 2

Station 27 WWSS 25 Canadian 15 WWSS 50 wwss 15 wwss 34 wwss 24 LRSM 55, mostly LRSM

Table 1

W. Longitude 116" 19' 1 9 17 56 24 30 22 12 28 49 27 21

Depth 666 m 819 1215 977 635 1161

Mean of magnitudes m. = 5.83 6.01 5.95 5.77 6.57 M. = 4.12 4.03 4.06

1200

71

Yield 70 kt 300 825 150

SD of magnitudes 0.47 0.40 0.39 0.41 0.54 0.32 0.33 0.22

52

Ulf Ericsson

site, all within 10 km from 37" 18' N, 116"23' W. The hypocentres were in rhyolite (Chartreuse, Halfbeak and Boxcar) or in tuff and below or very near the water table. The yields are measured in kiloton units of 4-18. 10l2J. They were radiochemically determined, as briefly described by Teller et al. (1968), and included prompt energy release due to gamma ray and neutron interaction with the hypocentre media. Higgins (1970) does not give the accuracy of the yields but quotes them as official for Chartreuse and as official preliminary and subject to future adjustment for the other five explosions. From the experimental yield ranges between f6 and f 15 per cent, determined by W. D. Richards for 12 explosions and reproduced by Teller et al. (1968), we can, however, conclude that the accuracy of such radiochemical yields is rather high when compared with most seismological energy measurements.

3. The magnitudes The magnitudes used for these explosions were the raw (uncorrected) magnitudes from 19 Canadian stations in the distance range from 12 to 49", as used by Basham (1969) in a paper on body and surface wave magnitudes from earthquakes and explosions. Magnitudes were lacking only in a few instances. The body wave magnitude m was calculated according to Gutenberg & Richter (1956) from vertical short period measurements on P,, or P , as explained by Basham (Zoc cit). The surface wave magnitude, here denoted by M , was calculated according to Vanek et al. (1962) from the maximum amplitude of the vertical movement of the Rayleigh waves. The data are reproduced in Table 3.

4. Themodel From inspection of the data a model simply linear in log,, W appeared as adequate for the functional relationship over the somewhat narrow log,,, W range from 1.8 to 3.1 and in a first step straight regression lines

Mij = a,'+&' log,, Wj+eij' mij = ai"+b,"log,, Wj+eij" were fitted for each station i and the explosions j measured there (The station and explosion indices i and j are defined in Tables 2 and 3). The ai' and a," ranged widely. The slopes bi' and bi" obtained did not show any correlation with distance or azimuth from source to station and neither did they appear to correlate with each other. According to their 95 per cent confidence ranges they were well compatible with the hypothesis of one slope b' and another slope b", common to all the 19 Canadian stations. Only a high b" for station GWC and a low b' for STJ were outside the 95 per cent range but well within 99 per cent. The station variances corresponding to the residuals eij' and eij"also appeared as uncorrelated with each other and according to Bartlett's test they were well compatible with the hypothesis of one eij"residual variance common to all stations and just compatible with the hypothesis of one eii' residual variance common to all stations. The offenders here were the somewhat large eij' residual variances at RES, SFA and OTT. These results implied three great simplifications of the m(W) and M ( W ) models for the Canadian stations, and possibly also for other stations: that the amplitude/ period ratios employed for magnitudes scale with source energy in essentially the same way at all stations but differently for short period P-wave movements and for vertical Rayleigh wave movements, that the corrections of station magnitudes to some common standard are independent of source energy and that the stochastic parts of the models for m ( W ) and M ( W ) possibly are nearly the same at all stations. The slopes b' and b" are the yield exponents in the scaling of amplitude/period

ALE MBC RES CMC FBC BLC YKC SCH GWC FSJ STJ FFC PHC HAL PNT SFA VIC SES O'IT

Station

1

5 - 5014. 54

--I-

6.0615.50/4.18 5.01/4*75

5.011-

5.5014.35 4-6814.64

5.411-

4.2114.85

-I-

5.011-

I

6

4 5

8 9 10 11 12 13 14 15 16 11 18 19

= 5.2214.70

=

5.0114.81 4.1014.75 4.1314.70 5.2514.53 5.1114.5 1 5.3814964

Explosion index m/M

2 3

Station index i=l 3

5.54/5.84

6.18/6.39 6.9815.92 6.50/6*01 5.8015.99 6.1616.09 6.5116.31

5.5515.41

5.6616 * 40 5.6716.60 5.1016.02 6.0115.96 5.6316.11 6.1416.04 6 * 5016.23 6.2216.11 5.8816.26 6.61j6.00 5.3915.91

5 -7816.12

6.08/5,81 6.83153 1 6.4515.70 6.2815.48 6.591541 6-1116.01

2 5.9115.61 5.5315.81 5.6616.31 5.4615.93 6*10/5,18 5.981980 6-2915.52 6.0715.81 5-18/5.71 5.6515.83 6.5315.53 5.9215.06

Table 3

5,5115.20 5.6215.44 5-2115.32 4.9615.62 6*11/5*31 5.1714.85 4.6015.15 5.5215.5 1 6.3615.02 6.0315.26 5-30/5*17 5.9115.21 5*10/5*60

5*58/5,41

49615.35 5.6015.21

5.0115.90

4 5.3915.22 5.3015 55 5

-15

46

5.5514-18

6-1114.59 5-3815.15 54914.50

-I-

4.8814.16 4.8715.16 4.5815.33 4.6314.15 5-1115.04 4.4815.00 5.3514.86 5.2514.81 4.8615.06 4.9015.09 5-5815.19 5.4614-59 4.8914.66

6.7116.06 -16.20 -16.27 6.1316.69 6.5516.31

-I-6.111-1-

6 6.2116.21 6.0416.31 6.2216.3 1 6.3116.28 6.0316,2315 *95/6.04 6.21/6-30 6-3016.39 6.5816.39 5.9916.22

w

u l

z

a

IF

%

5 e a

54

Ulf Ericsson

for yield. Their essential constancy between the stations is only a numerical convenience and not the consequence of some rigid geophysical law. This can be seen from the Fourier spectra of the experimental source functions of underground nuclear explosions, and their scaling law, e.g. in the papers by Carpenter et ul. (1962), by Werth & Herbst (1963) or by Haskell (1967). The spectra are cut off with a hump at high frequencies and increasing yields not only increase the amplitudes but also shift the spectrum and its cut off towards lower frequences. The overall result of an yield increase is an increase of low frequency amplitudes and a decrease of high frequency amplitudes. The ensuing changes in amplitude/period in the station record and thereby the numerical values of b' and b" depend on the detailed spectral shape of the Fourier transfer function from source function to record and this can well be different for the same wave phase recorded over different paths and with different instruments. In an attempt to obtain some b' and b" values also for other paths than from Nevada to Canada, the author searched the bulletins from stations UPP, KIR, NUR, KJN and KEV for m and M data from underground Nevada explosions with known yields. The b' and b" estimated from these data showed a tendency to be lower at these stations than in Canada, but the observations were too few for a significant numerical conclusion. The simplified functional relationships Mij = a,'+b' log,, Wj+eij' mij = ai"+b" log,, Wj+eij" with one slope b' and another slope b" were then adopted for the Canadian station. Refitting gave, with 95 per cent confidence ranges b' = 1*19f0-09 b" = 0.93f0.12.

The value b" = 0.93 is within the 95 per cent confidence level interval for all the 19 individual station slopes b,". The same holds for b' = 1.19 except for STJ, where it is only within the 99 per cent confidence range. The value b" = 0*93+0.12 is generally comparable with some of the slope values quoted in the articles referred to in the introduction and also with slopes obtained from m read off synthetic seismograms, computed by Kogeus at this institute for nuclear explosions in granite, with a wide range of yields and for instruments comparable to those in Canada. On b' nothing seems to have been published so far, but the value 1.19+0.09 compared well with other data shown by Evernden at the UMC 1970 conference in Stockholm on Geophysical Theory and Computers. From those data it also appeared that the slope holds down to yields below 10 kt. The means over the station intercepts ai' and ai" were a'. = 2.67 for M and u". = 3.49 for m. The differences 6.' = ai'-U' di" = a."-a" ranged widely, for M from -0.40 at FFC station to +0.35 at RES and for m from -0.49 at CMC to +0.81 at PNT. The numbers obtained are given in Table 4. The negative of these differences can be used as yield independent station corrections to a mean magnitude over the 19 stations. They are, as expected, well correlated with the Canadian station magnitude corrections obtained by Basham (1969). The and dimdo not show any obvious dependence on the azimuth from source to station. For ai"the correlation with the distance Ai from source to station is not significantly different from zero. The correlation coefficient between di' and log Ai is, however, nearly significantly different from zero and the regression of 6' on log Ai is 6i' = 0.49 log,, Ai-0.71

55

Yield dependence of magnitudes

Table 4 Station ALE MBC RES CMC FBC BLC YKC SCH GWC FSJ STJ FFC PHC HAL PNT SFA VIC SES OTT a,' = 2.67 a,"= 3.49

Station index i = l 2 3 4

6'' = u~'-u,' = -0.07

0.17 0-35 -0.01

-0.04 0.03 -0.13 0.07 0.07 0.16 0.02 -0.40 -0.21 0.29 -0.26 -0.09 -0.15

5

6 7 8 9 10 11 12 13 14 15 16 17 18 19 ~ ( 8 1 ' )= 0.19 ~(8,")= 0.37

s(e JI) = 0.19kO.03 s(el,")= 0 . 3 0 t 0 . 0 5

6," =

= -0.32

-0.30 -0.38 -0.49 -0.01 -0.20 0.07 0.09 -0.12 -0.34 0-36 0.04 -0.46 0.05

0.81 0.40 0.02

0.00

0.52

0.20

0.27

b' = 1'19+0.09 b" = 0.93k0.12

s3(e"1j)/s2(e'fj) = 2.34

where Aiis in degrees. A similar distance dependence can also be seen in the station corrections for M obtained by Basham (1969) from a larger data set, including the one used here and was likewise observed by Reagor et al. (1968) and by Jordan et al. (1970) with similar surface wave magnitudes from the Gasbuggy and Rulison explosions. It is a question of method whether the high magnitude values should be treated as regional anomalies or whether the magnitude definition for the region should be revised. Basham (1970b) has recently made such a revision. The differences Si' and Si" can also be seen as a stochastic sample of relative transfer functions from Pahute Mesa to various locations, the sampling having been done in the past by seismograph station site selection in a frozen field of stochastically distributed crust and mantle properties. The Si' and Si" appear to be uncorrelated with each other and according to the chi-square test at the 95 per cent level they show no significant deviation from normal distributions. Their estimated standard deviations are ~(6,')= 0.19 and ~(6,")= 0.37. The coefficients of correlation between the residuals eij'and eij"are important for joint uses of network measurements. In terms of an M and an m measurement at the same station being made by two different measurement channels, the matrix of correlation coefficients between the residuals in all channel pairs was calculated. Among the 703 non-diagonal coefficients only five were on the 99 per cent level significantly different from zero (e'(PNT)/e'(YKC), e'(PNT)/e'(SJT) and e"(OTT)/ e"(SFA) with positive coefficients and e'(GWC)/e"(STJ) and e'(PNT)/e"(SES) with negative coefficients), 31 channel pairs were undecided and 667 were not significantly different from zero at the 95 per cent level. The appearance of PNT in three possibly correlated pairs is curious, but in view of the rather few observations at each station the author regards the indicated correlations as spurious and takes all correlation coefficients to be essentially zero. The standard deviations common to all the residuals eij' and eijnat all stations

56

Ulf Ericsson

and their 95 per cent confidence ranges are s(eij’) = 0.19f0.03.

s(eij”)= 0.30f0.05.

The common variance of the surface wave magnitude residuals is clearly less than that for body wave magnitudes, as might have been expected from the wavelengths involved. One has s2(eij”)/s2(eij’) = 2.34. The standard deviations s(eij’)and s(eij”)measure the average variability, at constant yield, of the magnitude residuals when the individual ai’ and ai” or the differences di’ and din are known and applied. When this is not the case, the standard deviations ~ ( 6 ~and ’ ) s(di”) also contribute and the compounded standard deviations S’ = S”

J(0.19’+0*192) = 0.27

= J(O-3O2

+ 0.37’)

= 0.48

describe the variability of magnitude residuals at constant yield, if a Canadian station is picked at random and the station magnitude corrections are not known. The large value 0.48 explains something of the problems with unspecified P-wave magnitudes. The normal distribution describes the whole samples of eij’ and eij”only roughly. The elij pass the chi-square test for normalcy at the 95 per cent level, but the eij’ only at 98 per cent and both samples are disturbed by outliers, generating negative skewness and excess especially in eij”. The observations from each station appeared as too few for a meaningful determination of individual station variances, but for improved bulk representation, such as weighting over the network, the stations were somewhat arbitrarily split into two groups with low and high estimated station variances for eij’and two other similar groups for eij”.The groups and their estimated mean residual variances s12(eij’), shz (eij’), S12(eijn),

Shz(eij”)

etc. are given in Table 5. These groupings also make the split eij’ residuals and the low variance eij” residuals to pass the chi-square test for normalcy at 95 per cent and isolate the disturbances into the eij” high variance group, where the chi-square test is passed only at 98 per cent. The groupings of the stations do not follow any obvious geographical pattern. Stations PNT, SCH, GWC, MBC and HAL appear as having low variance in both kinds of magnitudes and stations STJ, CMC and BLC as having high variances in both kind of magnitudes. The residuals eij’ and eijnare physically generated by variations between records

Table 5 sl(e’lj)= 0.10

ALE, MBC, YKC, SCH, GWC, FSJ, FFC, PHC, HAL, PNT, VIC

sl(e’,,) = 0.11

,,

RES, CMC, FBC, BLC, STJ, SFA, SES, OTT

sb2(e’,J= 0.067 sh2(dij) = 0.036

MBC, RES, FBC, SCH, GWC, HAL, PNT, SFA, SES, OTT

s12(e”ij)= 0.033

ALE, CMC, BLC, YKC, FSJ, STJ, FFC, PHC, VIC

No signif. deviation from normalcy

S,’(&”,,)

= 0.025

sh(e’)= 0.26 = 0.19

sh(&’fj)

s1(e”,,) = 0.18 S I ( E ” ~ J )= 0.16

s,(e”,,) = 0.37 sh(&”iJ) = 0.30

,, ,, ,, I,

Nearly signif. dev. from normalcy

57

Yield dependence of magnitudes

Table 6 1 -0.24 0.20

2 0.07 0.14

3 -0.01 0.09

4 0.08 0.09

0.07 0.17

6 -0.06 0.16

e,J” = -0.09 SJ(E~J”) = 0.26

0.26 0-13

-0.14 0.26

0.09 0.35

-0.04 0.38

-0.02 0.37

j = e.,‘ = SJ(E1;)

=

=

0 . 11

&I)’

= eiJ’-f?.J’

s(e.J”) =

0.13

Eij”

ell“= -e.J”

s(e.J?

5

s(&IJ’)= 0.16 S(&]t) = 0.27

in the data analysis, by earth noise and by variations between explosions of the energy coupling in the source area. They also contain such higher order terms as might have been neglected by the simple linear model employed here. The energy coupling will depend on the details of the hypocentre geology and on how small and deep down the emplacement cavity is. The present data are rather few per station and do not offer good opportunities to separate such details. We can, however, obtain measures of the bias in coupling of each explosion by taking the averages e.j’ and e.j” over stations of the residuals for each explosion. These averages are given in Table 6 and show how the Chartreuse explosion made generally low M and Halfbeak generally high m,for the other cases the coupling of the explosions did not deviate much from the average. With the e.j‘ and e.j” and the unbiased residuals E..’ = eij’-e.j’ qj”= e..”-e IJ

.”

.J

one can then separate the variances sz(eij’) and s2(eij”)into variances s2(e.j’) and s2(e.j”) between explosions and variances s ’ ( E , ~ ’ ) and sz(cij”)between stations as s2(eij‘)= s2(e.j’)+s2(eij’)

s2(eij”)= sZ(e.j”)+s2(Eij”).

The present data give (0.19)’ = (0*11)’+(0-16)’

(0.30)’ = (0.13)’+ (0.27)’. The variation between explosions was less than between all stations. The s(cij’) and s(Eij”)obtained after the earlier discussed split into groups of high variance and low variance stations are given in Table 5. Both the unsplit and split cij’ and cijX passed the 95 per cent level chi-square test for normalcy, except the E ~ from ~ ” the high variance stations. The standard deviations sj(cij’) and sj(cijn) of the unbiased residuals for each explosion measure the variabilities of the noise and radiation pattern in the individual explosions, as obtained in the azimuth quadrant in Nevada occupied by the Canadian stations. The estimated values are given in Table 6, the highest and lowest sj of both kinds appear to be signjficantly different on the 95 per cent level and again the m magnitudes show more variability than the M magnitudes. The present data do not provide for separation of radiation pattern from analysis and earth noise.

5. Magnitudes expected from given yields The magnitudes expected, for a given yield, at any of the 19 stations above are directly obtained from the functional relationships in our models for the station in question. One application of such expectations could be to a proposed ban on underground nuclear explosions producing more than some agreed threshold magnitude

58

Ulf Ericsson

One intention with the threshold level would be to make the underground events with magnitudes above the threshold strong enough to be easily identifiable as nuclear explosions or earthquakes. From the foregoing it should be clear that the definition of such a magnitude threshold should involve careful definitions of the kinds of magnitudes, the stations and the station corrections to be used, possibly with different stations for different source areas, of how to compound measurements from the stations and how to deal with stations that fail to record. The matter would be quite involved and even so it appears from the stochastic part of our model, that there will always be a non-zero probability to overstep the threshold by an explosion expected to be well below it. In addition there will be some overstepping by above threshold earthquakes misidentified as explosions. A threshold ban would hence also need provisions for how to deal with really and apparently overstepping events. An explosions-making agency will calculate its probability for violation by first using some model fo r the probability of having an explosion of given yield detected at given stations, some model like ours to obtain the probability of then overstepping the threshold and finally some model for the probability of having the explosion identified as an explosion. Decisions to explode would be according to some politically determined and probably very low probability for overstepping. A distinct transition of probabilities over the threshold would minimize occasions for technical and political controversy. This would be aided if overstepping probabilities were compounded from several stations, if surface wave magnitudes M with their inherently small residual variances were used as much as possible and if generally low variance stations were preferred before high variance stations. For a simple example, involving only the probabilities for overstepping a threshold, we take the politically much discussed magnitude threshold m = 4.75 and apply it first to the MBC station. According to our model, now somewhat extrapolated to low yields, this threshold means that 50 per cent of 48 kt explosions in Pahute Mesa are expected to give m below 4.75 at MBC and 50 per cent above. The corresponding 50 per cent surface wave magnitude at MBC is M = 4.84. We then assign the standard deviation s(dij) = 0-1and s(e'lij) = 0.18 to the MBC station and obtain the following probabilities P for overstepping the MBC station magnitude thresholds m = 4-75 and M = 4.84:

P = 50% 30 10 5

W' = 48 kt

1

43 37 35 30

0.1

26

W " = 48 kt 38 27 23 17 12

This exemplifies the difference in transition range for the yield. In practice, however, the important considerations would not be made at the 50 per cent level but at some level more like 0.1 per cent and we see that in our example above the yield thresholds are different at that level. Such rather undesirable situations are avoided if the threshold is redefined altogether to be defined in kilotons at the source with an associated low overstepping probability equal at all stations in the network. An example would be 20 kt in Pahute Mesa with 0.1 per cent overstepping probability. If the low level and high level variances for the e residuals in Table 5 are adopted and the extrapolation of our model is accepted for the exercise, the corresponding magnitude thresholds would be as given in Table 7. In a refined calculation of such thresholds the calibration quality at each station would be taken into account. The simple exercise here only illustrates one practical application of our model for network m ( W ) and M ( W ) .

59

Yield dependence of magnitudes

Table 7

Station ALE MBC RES CMC FBC BLC YKC SCH GWC FSJ STJ FFC PHC HAL PNT SFA VIC SES OTT

0.1 per cent threshold magnitudes for 20 kt in Pahute Mesa in Nevada M = 4.47 m = 5.27 4.71 5.37 5.01 4.99 5-05 4.41 4.62 4.62 4.70 5.04 4.62 4.33 4.83 4.28 4.93 4.39 5.03 5.22

4.71 4.62 5.09 5.00 5.39 5.66 5.09 4.89 5.25 5.94 5.62 5.12 5.05 5.81 5.40 5.60 6.10 5.28

6. Estimates of yields from magnitudes

One interesting technical application of our model is to the estimation of explosion yields from measured magnitudes. The normally distributed residuals in our model immediately make available the regular statistical formulae for maximum likelihood estimates and log,, @ of the log yield K = log,, W , with their associated confidence ranges. They are weighted means over the successfully recording stations or channels of the network. We will give some of the expressions in a form that always extends over the whole network but with provision for failing channels. The maximum likelihood estimates of K are unbiased and have minimum variance. They of course minimize only the stochastic variability and assume that the conditions under calibration hold also under application. They do not eliminate systematic errors due to, say, exceptionally porous hypocentre media. The yield estimates # are not true yields, only 'seismic' or 'equivalent' yields, the latter referring to source conditions equivalent to those during calibration. A derivation of the maximum likelihood estimate in the case of one channel with a model like ours, is found in Graybill's (1961) treatise on statistical models. For the multichannel case and with our notation one obtains Iog PP = R

Ri,'

= (Mio- u.' -6'i)/bi'

Kion= (mi, - u". -ai")/bi"

where the Mi, and mi,are the magnitude observations obtained from the unknown yield. The success parameters and f i 0 are unity or zero, depending on whether the channel in question succeeded to provide an observation or not, and all sums extend over all stations in the network. The a.' u." di' bi" bi' bi" si(eij') and si(eij") +Iio

60

Ulf Ericsson

are obtained from earlier calibrations on known yields. For the standard deviations, however, only the ratios are required. With b' being larger than b" and s(eij') being smaller than s(eij"), as in Table 4, the weight of M-magnitudesis about four times larger than the weight of rn-magnitudes in the estimation of yield. For yield estimates l?' or I@" from M or m measurements only, the expressions simplify accordingly. Other simplifications occur if the slopes b' and b" or the estimated variances are equal in groups, as appears to be the case in our present material. The locus of a fixed maximum likelihood estimate of I? is a hyperplane in the Mi mi space. In the common bichannel case, as found in many observatories with one vertical short period instrument and one vertical long period instrument, this reduces to a linear relation between only two variables. With the indices dropped it is

A family of such straight lines in an m(M)diagram provides a convenient nomogram for optimum bichannel estimates of yield. Fig. 1 shows such a W ( M ,m) nomogram for Pahute Mesa explosions measured at the MBC station. The constants for the functional relationships were taken from Table 4, the variances were assumed to be as given in Table 5. In the case of non-zero correlation p between the bichannel residuals e' and e" the squared fractions in the expression above for log l? would have to be changed

x

PAHUTE/MBC

m

7 -

1000 k t

6 -

5 -

4-

3-

24 2

I

I

I

I

3

4

5

6

I

I

7

8

M

FIG.1. Maximum likelihood estimates fl of the yield of nuclear explosions in Pahute Mesa, obtained from m and M measured at MBC station.

Yield dependence of magnitudes

61

into

and 2

b’

b”

(&)-Ps(e’).s(e”)*

In a multichannel case with non-zero correlation between channel residuals, the maximum likelihood estimate R is conveniently written in matrix notation: hV}R = (W- W T

UJ}

(b}T{COV}-~ (b}

where {m}

is the observed column vector of magnitudes mi,

@I

is the column vector of ai, obtained by Calibration

(4

the same for the slopes

{cov}-’

is the inverse of the matrix {cov} of covariances between residuals, known beforehand or obtained by calibration.

and

and where the indices i now indicate channels. The explicit 2-formulae above would reappear, if one would use the covariance matrix to transform from the original coordinates Mi m ito a proper orthogonal set of uncorrelated variables. Such orthogonal ‘ magnitudes ’ would have somewhat more physical meaning than the current ad hoc defined magnitudes. The present material appears, however, to have esentially uncorrelated residuals, so there is no immediate reason to enlarge on the use of such transformations. Confidence limits at the ct level for K can be obtained as the K-roots of the expression

where I is the number of channels giving at each channel the logyield estimates for the unknown logyield, ai and b, are the channel intercepts and slopes obtained by calibration with J events with known logyields Kij and Rij = (mij-ui)/bi the associated channel estimates. F , is Snedecor’s F at the a level with degrees of freedom I and

Ri, = (mio-ai) b,-l

The Ki. are the means of Kij in a channel, the 4ijare the success parameters in the calibration and the ci relative values for the standard deviations of residuals in each channel. The ci are supposed to be population values or very good estimates of them

62

Ulf Ericsson

and so make these confidence ranges somewhat optimistic. In the case with just one channel, Z = 1, the expression reduces to (Graybill 1961):

K = R , i f , , , . s . l+J-’+(K-K.)‘

)-l)+.

j=I

This expression shows, with dropped station indices i, how the confidence range is influenced by the estimated standard deviation s of the residuals, the number J of calibrating explosions, the relative position K - K . of the unknown explosion K and the mean calibration explosion K. and by the parameter X(Kj-K.)’ of the calibration explosions. The distribution is now the two tailed Student’s t on the M level and with 5-2 degrees of freedom. As an example, Fig. 2 shows the 95 per cent confidence ranges for Pahute Mesa logyields estimated from the MBC station M on the basis of the MBC calibration data in Tables 2, 3, 4 and 5.

7. The use of m,M observations for event identification Simultaneous measurements of m and M have proved to be quite effective for the identification of earthquakes and underground nuclear explosions, a task of political interest in a ban on underground nuclear explosions. At constant M the earthquakes mostly produce smaller m than the nuclear explosions. The results on the detailed behaviour of m,M have, however, been rather conflicting, see the papers by Marshall et al. (1966), Capon, Greenfield & Lacoss (1967), Thirlaway (1968), the SIPRI (1968) report, Basham (1969), Lieberman & Pomeroy (1969) and Pasechnik et al. (1970). One important point is if and how the distributions of explosion and earthquake m, M points converge and degrade the identification when events get weaker. The author developed a method (Ericsson 1968, 1970) for assessing the political usefulness of particular identification methods and for that method a description of the statistical behaviour of identification parameters is essential. Our model provides such a description for m(M) observations on underground nuclear explosions.

/

1..1-11

4 -

3 -

’

*’

/’ J’

2 ‘.

0

’

I

/

/ /

+ I

2

J

/

3

4

I

I

5

6

I t -

7

8

ni

FIO.2. Maximum likelihood estimate k,, and 95 per cent confidence intervals of log W for nuclear explosions in Pahute Mesa, obtained from M measured at MBC station.

Yield dependenee of magnitudes

63

If we consider a series of explosions all with equal yield W, our model says that the my M obtained at a bichannel observatory i will be distributed according to a binormal distribution in the m, M plane. The ellipses of equal probability or likelihood will be tilted, if there is correlation between the residuals eij’ and eijN. From our present material, however, we found that the correlation was essentially zero, so that these ellipses will be parallel with the m and M axes. They will be centred at M

= ai‘+b’

m

=

log W

ailt+b” log W

and for other yields they will be centred on the straight line

m

+ (b”/b‘)(M- uil).

= air’

The m residuals d” from this line of mean or expected m, M or functional relationship between m and M will be roughly normally distributed with variance

~’(drr)

=

+

sz(eijl’) s’(eij’)(b”/b’)’,

when the correlation is zero. With the earlier estimated slopes we obtain b“/b‘ = 0-78 with the standard deviation 0-06. From the variances in Table 4 we obtain ~(d”)’ = (0.30)2. This means that an explosion generated m, M at a Canadian station has a chance of about 10 per cent to be more than 0.42 units below the mean m(M) line for that station. Whether this is sufficient or not for identification depends on where the earthquakes are distributed in the m(M) diagram, how frequent they are and what the political requirements on the probability for correct identification of explosions and on the avoidance of false alarms are. This has been discussed in eartier publications by the author (Ericsson 1968, 1970). A small s(d”) of course aids discrimination. One particular conclusion from our model for m and M from the Canadian stations is that the position of the m(M) line will be rather different for different m, M channel conibinations. The line can shift as much as 1.9 units in the m-direction, between the extreme cases with m from CMC and M from RES, being 0.76 units below the mean Pahute-Canadian Iine m = 0.78 M + 1.41 and nz from PNT and M from FFC, being 1.12 units above the mean line. Mixing of uncorrected my M from various stations adds to the variance s2(d”) and degrades ~) the discrimination capability. Compounding with the earlier determined ~ ( 6 ’ and ~ ( 6 ’ ’ shows ~) that s(d”) would increase from 0.3 to 0.5 in such a case of channel mixing. The effect is avoided if station corrections are applied before mixing, in our notation one would plot mi-&“ against M i -ai’ and would obtain a scattergram around the line above. The variance of station or bichannel residuals would then be sz(d”). Basham (1969) plotted m, M obtained as means, over the Canadian network, of station corrected magnitudes for many Nevada explosions. This of course corresponds nearly to the procedure suggested here and the taking of means reduces the residual variance below the bichannel value sz(d ). Such residual variances of means would be even further minimized and the discrimination sharpness would be increased if the network means were taken as weighted means according to the grouping and variances given in Table 5. Basham (1969) calculated mean magnitudes for the same Pahute Mesa explosions as studied in this paper and also for some other explosions at the Nevada Test Site. For the Pahute Mesa explosions he used the same station magnitudes as used here. In addition his station magnitude corrections were not very different from the -hi’ and -di“ obtained here, so that his unweighted means mCANand McAN over the Canadian network can be compared with means from our model. Basham’s mCAN and M C A N are reproduced in Fig. 3, with explosions in Pahute Mesa and under

64

Ulf Ericsson

Yucca Flat differently marked. Our line, when applied to the mCAN and MCAN, fits the Pahute Mesa explosions quite well and also the mCAN, MCAN point observed by Basham, Weichert & Anglin (1970) for the explosion Benham in Pahute Mesa. Benham was not used in the determination of our model and thus serves as a check. We also see that the Yucca Flat explosions generally fall 0-3 units above our Pahute Mesa line and nearly parallel to it. This shows how the u.’ and u.’’ of our model are more sensitive to the location of the source region than the ratio b“/b’. Our line for Pahute Mesa and a parallel line through the Yucca Flat points differ markedly from 7 by Basham (1969) to cover both the line MC-N = 1.43 m c ~ ~ - 2 * 8determined source areas. This is important for the identification of weak events. It is easy to see that the linear function D = m - (b”/b’)M + (b”/b’)a‘ - a’’ of m and M is independent of K or W,as long as the functional relationship m(M) is linear. Our model thus provides us with a natural candidate D for a discriminant with the desirable property of being independent of the yield W. In most cases of practical application the functional relationship m = a”+ (b”/b’)(M-u f ) and the variances would have to be estimated from the m yM points without the help of known explosion yields. Such is always the case with earthquakes, where the yields or energies appear quite unaccessible. The estimation of the functional relationship and the variances in such cases is the subject of a paper in preparation. 8. Conclusions

The yield dependence, at each of a network of 19 Canadian seismograph stations, of the body wave magnitudes and surface wave magnitudes generated by six nuclear explosions in the Pahute Mesa of Nevada Test Site and with independently meas1 red

mean

Canadian means

I

6 -

5 -

/

4 -

2

0

3

4

5

6

7

8

Mcan

FIG. 3. Functional relationship mCAN(MCAN) for nuclear explosions in Pahute Mesa, marked by filled circles and by the cross (explosion Benham). Open circles are nuclear explosions in alluvium or tuff under Yucca Flat. Magnitudes mCAN and M,,, are Canadian means. They and the dashed line were obtained by Basham (1969, 1970a).

Yield dependence of magnitudes

65

yields from 70 to 1200kt, were described by functional or deterministic relationships simply linear in the logyield plus stochastic parts independent of yield. In the functional parts of the model the intercepts varied widely between stations, by more than one unit for the body wave magnitudes. The intercepts were nearly normally distributed, with standard deviations of about 0.2 and 0.4 for surface and body wave magnitude intercepts respectively. They did not show any obvious correlation with each other or with geography, except a nearly significant positive correlation of the surface wave magnitude intercepts with the distance from Nevada Test Site to the station. The slopes, also being the exponents for yield scaling of amplitude/period, were essentially equal at all stations but larger for surface wave magnitudes than for body wave magnitudes. Their 95 per cent confidence ranges were found to be 1.10 to 1.28 for surface waves and 0.81 to 1.05 for body waves, with point estimates 1.19 and 0.93 respectively. These values for body wave magnitudes can be compared with slopes in other publications but the surface wave slope or yield scaling exponent seems to be new. In the stochastic part of the model, the magnitude residuals appeared as uncorrelated between stations. Their standard deviations at all the stations can be described by one value, 0.2, for surface wave magnitude residuals and by a significantly larger value, 0.3, for body wave magnitude residuals. The common distribution of these residuals are then only roughly normal. The representation is, however, improved if the stations are split into certain groups with low and high residual variances. The distributions of the residuals than become nearly normal. The model gives access to separate estimates of the variation between explosions and of the variation between stations of magnitude generation and ' energy ' coupling to earth. The latter variation appeared to be the larger. The constancy of the station slopes implies that the station magnitude corrections to some common standard are independent of explosion strength. The obtained corrections of the station magnitudes to a mean over these Canadian stations, compared well with station corrections obtained by Basham (1969) from these and other nuclear explosions in Nevada, New Mexico and Colorado. If, for a given yield, the magnitude at a randomly picked Canadian station is estimated without knowledge about the intercept for that station, the standard deviations of the intercepts and of the residuals compound and produce overall standard deviations of about 0.3 and 0-5 for surface and body wave magnitudes respectively. This explains the often observed large magnitude differences between stations generally. In regard to applications, the wide range of station intercepts and the unequal residual variances for body and surface wave magnitudes imply that a threshold for a ban on underground nuclear testing should not be defined by a single magnitude but by a yield and by associated magnitude thresholds individual to stations and defined for some low and common overstepping probability. For application to the estimation of unknown yields from the magnitudes measured in a network, the maximum likelihood estimate by a properly weighted mean becomes available and also a confidence range that takes the quality of network calibration into account. Because of the lower standard deviation of their residuals, the surface wave magnitudes are more precise for yield estimation. In application to the identification of explosions and earthquakes by simultaneous measurements of body and surface wave magnitudes, the model explains how the position, slope and variability of such measurements on explosions are generated. For individual Canadian stations the positions of the linear functional relationships between body and surface wave magnitudes vary widely, but they remain parallel. For a station with known intercepts the standard deviation of the m residuals against its functional relationship is about 0.3 and larger when stations are mixed without E

66

Ulf Ericsson

station correction. The mean of station corrected functional relationships compares well with unweighted but station corrected Canadian mean magnitudes for explosions in Pahute Mesa but shows a difference for explosions under the Yucca Flat in the Nevada Test Site, presumably because of shifted intercepts. The model also immediately provides expressionslinear in body and surface wave magnitudes but independent of the yield. These should be investigated for use as yield independent discriminants between explosions and earthquakes.

Acknowledgments The author has been much aided by the programming and computational work by K. Kogeus, G. Hoernstroem and I. Eriksson, by the secretarial work of G . Esko and by discussions with 0. Dahlman and H. Israelson at this institute. The investigation was made possible by the extensive magnitude data kindly supplied by Dr P. Basham at the Earth Physics Branch in Ottawa, Ontario, of the Canadian Department of Energy, Mines and Resources.

Research Institute of National Defence Stockholm, Sweden

Notations c(

confidence level.

ai' and ail' estimated intercepts (and decadic logarithms of coupling factors) at station i. a,' and a," means of station intercepts ai' and ai". a and a, intercepts of channel i.

bi' and bi" estimated slopes (and yield scaling exponents of amplitude/period) at station i.

b, b' and b" estimated slopes for channels generally and for M and m channels, respectively. Ci

relative value of the standard deviation of magnitude residuals in channel i.

cov co-variance between magnitude residuals. di' and di" a,' -a.' and afIi-a".. d

I'

residual of m from functional relationship m(M).

D discriminant. hi geocentric distance from Nevada Test Site to station i, in degrees. eij' and eij" magnitude residuals against the functional relationships M ( W )and m ( W ) at station i for explosion j. e,j' and e.j'' means over stations of eij' and eij". eij' and cij" unbiased residuals eij' -e.j' and eijf'-edj". F,

Snedecor's F at confidence level a.

Yield dependence of magnitudes 67 success parameters, equal to 1 or 0 depending on whether the measuring station or channel indicated by superscript and index provided a magnitude or not.

I

number of channels.

J

number of calibration explosions.

K

the logyield log,, W.

Kij logyield of calibration explosions j for channel i. Ki. mean, over calibration explosions, of Kij.

Izij

estimate of Kij from mij.

M or Mij surface wave magnitude measured from the maximum amplitude of of the vertical Rayleigh surface wave movement at station i due to explosion j.

m or mij body wave magnitude from vertical short period P,, or P wave movement, at station i due to explosion j. In paragraph 6 the mij refer to body or surface wave magnitude as appropriate to measurement channel i. unweighted means over Canadian stations of station corrected m and M . compressional body wave running along the interface of the Earth's crust and mantle. P

compressional body wave running through the Earth's mantle.

P probability. coefficient of conelation. estimated standard deviation. estimated standard deviations of quantity between parentheses. standard deviations of M and m residuals estimated by use of mean intercepts a'. and a".. estimated standard deviations for station groups with low (1) and high (h) magnitude residual variances. sj(&ij') estimated standard deviations of the unbiased residuals for explosion and sj(qj") j.

W

explosion yield or energy, in kiloton units of 4.18. 10l2joule each.

@ @' and @" maximum likelihood estimates of the yield from all measurements jointly, from M measurements only and from m measurements only. these brackets denote a matrix, its transpose and its inverse. the cap on a quantity, like 8, denotes a maximum likelihood estimate. superscripts ' and " refer to surface wave and body wave magnitudes respectively.

68

Ulf Ericsson

o index o (single or second) indicates an explosion with unknown yield W, and logyield KO,to be estimated. j

index j or the second index of a double index indicates the explosion j, in the present application specified in Table 2.

i

index i or the first index in a double index indicates a seismograph station, as specified in Table 3. In paragraph 6, however, index i indicates a channel measuring m or M , as the case may be.

.

the dot index . indicates the arithmetic mean over the index it replaces.

References Basham, P., 1969. Geophys. J. R. astr. SOC., 17, 1-13. Basham, P., 1970a. Can. J . Earth Sci., 7, 531-534. Basham, P., 1970b. Contribution from the Earth Physics Branch, Department of Energy, Ottawa, Ontario. Basham, P. W., Weichert, D. H. & Anglin, F. M., 1970. J . geophys. Res., 75, 15451556. Capon, J., Greenfield, R. & Lacoss, R., 1967. Technical Note 1967-50, Lincoln Laboratory, Lexington, Mass. Carpenter, E. W., Savill, R. A. & Wright, J. K., 1962. Geophys J. R. astr. SOC.,6, 426-440. Ericsson, U., 1968. Report C 4338-20, Research Institute of National Defence, Stockholm. Ericsson, U., 1970. Bull. seism. SOC.Am., 60, 1521. Evernden, J., 1969. Trans. Am. geophys. Un., 50, 247. Evernden, J., 1970. J . geophys. Res., 75, 1028-1032. Graybill, F., 1961. An introduction to linear statistical models, McGraw-Hill, New York. Green, P., Frosch, R. & Romney, C., 1965. Proc. I.E.E.E., 53, 1822. Gutenberg, B. & Richter, C. F., 1956. Ann. Geofis., 9, 1 . Hankins, D., 1964, Proceedings of the Third Plowshare Symposium, TID-7695, 195206. Haskell, N., 1967. J . geophys. Res., 72, 2583-2587. Am., 59, 2283-2293. Hays, W., 1969. Bull. seism. SOC. Higgins, G., 1970. Report UCRL-72346, Lawrence Radiation Laboratory, Livermore, California. Jensen, O., Ellis, R. & Russell, R., 1966. Earth Planet. Sci. Letts., 1, 211-221. Jordan, J., Black, R. & Bates, C., 1965. Bull. seism. SOC.Am., 55, 693-720. Jordan, J., 1970. U.S. Coast and Geodetic Survey report Seismic data from Rulison. Lambert, D. G.,von Seggern, D. H., Alexander, S . S . & Galat, G. A., 1969. Seismic Data Laboratory Report 234, Alexandria, Virginia. Latter, A. L. & Martinelli, E. A., 1959. Phys. Fluids, 2, 280-282. Lieberman, R. & Pomeroy, P., 1969. J . geophys. Res., 79, 1575-1590. Marshall, P. D., Carpenter, E., Douglas, A. & Young, J., 1966. AWRE-report 0-67/66, Aldermaston, Berkshire. McGarr, A., 1969. Bull. seism. SOC.SOC., Am., 59, 1307-1334. Pasechnik, I. P., Kogan, S. D., Sultanow, D. D. & Tsibulskii, V. I., 1960. Trans. 0. Schmidt Znst. Geophysics, 15(182), Moscow. Pasechnik, I., Dashkov, G., Polikarpova, I. & Gamburtseva, N. G., 1970. Earth Physics, 1, 28-36.

Yield dependence of magnitudes

69

Reagor, B. G., Gordon, D. W. & Jordan, J. N., 1968. ESSA report PNE-G-32, Rockville, Maryland. Riznichenko, Y., 1960. Trans. 0. Schmidt Institute of Geophysics, 15, (182), Moscow. Romney, C., 1959. J. geopliys. Res., 64, 1489-1498. Romney, C., 1962. VESIAC Special Report 4410-36-X, University of Michigan, Ann Arbor, Michigan. Romney, C., 1963. Transcript of hearings of the Joint Committee on Atomic Energy, 88th U S Congress, March 5-12, 1963. SIPRI report on Seismic Methods for Monitoring Underground Explosions, 1968. Stockhclm International Peace Research Institute. Springer, D., 1966. Bull. seism. SOC.Am., 56, 861-876. Teller, E., Talley, W., Higgins, G. & Johnson, G., 1968. The constructive uses of nuclear explosives, McGraw-Hill, New York. Thirlaway, H. I. S., 1968. Contemp. Physics, 9, 17-31. Werth, G. C., Herbst, R. F. & Springer, D. L., 1962. J. geophys. Res., 67, 1587-1610. Werth, G. C. & Herbst, R. F., 1963. J. geophys. Res., 68, 1463-1475. Vanek, J., Zatopek, A., Karnik, V., Kondorskaya, N. V., Riznichenko, Y . Y., Savrensky, E. F., Solovjev, S. L. & Shebalin, N. V., 1962. Bull. (Izvest.) Acud. Sci. U.S.S.R., Geophys. Ser., 2, 108.