Notes on moment tensor inversion programs INVRAD.FOR INVDC.FOR SYNRAD.FOR by John E. Ebel, Dec., 1988. The above three programs were written in the summer of 1987 at the Geophysical Institute of the University of Karlsruhe in West Germany to test the idea that the amplitudes of the direct arrivals from small earthquakes could be inverted for the source focal mechanism. Each program performs a different but related reduction of the input amplitude data, as follows: INVRAD.FOR - Inverts the direct P, SV and SH amplitudes to find a traceless moment tensor using the generalized inverse method as outlined in Aki & Richards. The program finds the largest double-couple component of the resulting traceless moment tensor by subtracting out a CLVD component and then returns focal mechanism, seismic moment and error information. Since the amplitudes are directly proportional to the moment tensor components, the inversion is non-iterative and is a least squares solution if 5 or more amplitude readings are entered to solve for the 5 independent moment-tensor components. INVDC.FOR - Inverts the direct P, SV and SH amplitudes to find the best constrained double-couple solution for the data. The generalized inverse method as outlined in Aki and Richards is used. Since the double-couple formulation for the amplitudes is non-linear, the program needs a starting focal mechanism (strike, dip and rake) and then attempts to iteratively find the best-fitting focal mechanism solution based on the data and the starting solution. The best focal mechanism, seismic moment and error information are listed. SYNRAD.FOR - Solves the forward problem of calculating the amplitudes of direct P, SV and SH arrivals at selected stations given a focal mechanism and seismic moment. Amplitude data are entered, and the fit of the predicted amplitudes to the observed amplitudes are calculated. In general, I use INVRAD to do initial inversions on the data, INVDC to see if a somewhat better constrained double-couple solution can be found using the solution from INVRAD as the starting solution, and SYNRAD to test any additional poorer quality data which I may have. Since INVRAD solves for a moment tensor and then extracts a double-couple solution from that moment tensor, the fit of the observed amplitudes to the predicted amplitudes from double-couple solution is not as good as the fit of the observed amplitudes to predicted amplitudes from least-squares moment tensor solution. For events with large CLVD components, differences in the amplitudes between the observations and those predicted by the double-couple solution can be quite large for individual data points. The program INVDC was written to try to improve the fits from INVRAD by constraining the solution to be a double couple, although in my limited experience to date, I have found the improvements to be quite minor. The program SYNRAD was written to test individual data points of questionable quality using solutions found from INVRAD and INVDC. Since INVRAD and INVDC give least-squares solutions for 5 or more data points, one or two bad data points can badly throw off the solution. My practice so far, again with limited experience, has been to use only amplitude readings of high confidence in initial runs using INVRAD and INVDC, then to use the solutions from these programs in SYNRAD to test my less confident amplitude readings. All of the solutions in my paper on the Black Forest earthquakes were found using INVRAD. Program assumptions A number of assumptions and simplifications are built into the present versions of these programs. First, the programs, as written, can only handle direct P and S waves (upgoing) from the earthquake source. All rays must intercept all velocity discontinuities at pre-critical angles. The subroutines which calculate the ray path, such as PATHDEF and NEWTTME, were initially programed with the idea in mind of allowing more complicated ray paths (such as reflections, phase conversions, etc.), but that has not been completed in this version of these programs. Second, the programs as written correct the ray amplitudes only for 1/R geometric spreading and not for anelastic attenuation or for the proper transmission coefficients at velocity interfaces. The first assumption should be acceptible in source areas of high Q. If reflected phases or converted phases are to be programed in, then the proper reflection/transmission coefficients must be included. Third, the crustal model must have planar layers, but it is allowed to have vertical velocity gradients within any or all of the layers. The velocity model resides in its own file where each line gives the depth (in km) and velocity (in km/sec) for the top or bottom of a crustal layer, starting with the values at the free surface. The format for each line is 2F10.4. A sample velocity model is: 0.00 5.0 2.00 6.0 2.00 6.0 10.0 6.0 10.0 6.6 20.0 6.9 20.0 7.1 34.0 7.2 34.0 8.1 This model has gradients in the top layer (5.0 km./sec to 6.0 km/sec over 2 km depth), then a constant velocity (6.0 km/sec) layer from 2 to 10 km. At 10 km is a velocity discontinuity, with a jump from 6.0 km/sec to 6.6 km/sec. Between 10 km and 20 km there is a velocity gradient from 6.6 km/sec to 6.9 km/sec. At 20 km, the velocity jumps to 7.1 km/sec and from there to 34 km there is a slight velocity gradient. Finally, at 34 km the velocity jumps from 7.2 km/sec to 8.1 km/sec. The S-wave velocity structure is assumed to be the same as the P-wave structure except for the usual factor of 1.732. Fourth, the parameter TAREA is the area of the farfield time function, and this must either be known or assumed in order to calculate correctly scaled theoretical amplitudes (TAREA scales the source time function to unit area). I have set it to .01 sec, which assumes essentially a delta function as the source time function shape. This is probably a pretty good value for small earthquakes (less than about 2.2 or so), but should be larger for larger events. Changing the value of TAREA by some factor C will change the theoretical (predicted) amplitudes for a given seismic moment computed in the programs by 1/C. Program Notes 1. The number of amplitude readings is presently limited to 50 by the dimensions in the main programs. Also, the maximum number of layers which can be traversed by a ray is limited to 50 by the dimensions in the main programs. 2. The programs were originally written on an HP-1000 computer where the screen read and write units were 1. These versions were modified for a VAX VMS machine where the read and write units are 5. You can modify these if need be. Also, the HP computer had a special statement EMA to handle larger programs and arrays. The EMA statements themselves have been deleted from the programs, but some of the statements and comments still refer to the coding necessary for the EMA format on the HP-1000. 3. The units for variables in the program are as follows: all amplitudes in microns (including estimates of the data errors needed for the variance calculations), velocities in km/sec, density in g/cm**3, strike, dip and rake in degrees, moment in dyne-cm (generally x 10**16), azimuths in degrees and all distances in km. The density of the source region which is set in the beginning of the program is 2.7 g/cm**3. 4. In the program INVDC the user is asked for a multiplicative factor for the model steps. This factor, which should be set between 0 and 1, is multiplied into all of the calculated perturbations to the model parameters as calculated in the inversion subroutine. The purpose of the factor is to control the inversion to be sure it does not overshoot a local minimum into some other, undesirable part of the solution space. 5. The programs contain extra subroutines which are not called or needed in one or more of the individual programs. I include them partly out of laziness and partly because they contain comments which may help explain how the subroutines work. 6. The programs do not check for rays past critical angle except to determine if SV rays are past critical angle at the free surface. If so, the programs list those rays past critical angle, and then they bomb. At this point the best remedy is to delete the offending station readings from the data list and rerun the programs. Note: The programs were revised slightly since December, 1988 when an initial distribution of the programs was made. Several bugs were found and corrected, and some of these bugs could give spurious results under some particular circumstances. I cannot guarantee that all errors are out of these programs, and I would appreciate it if you would let me know about any other bugs which you may find. I would also like to know of data sets where the programs were found to work well, and those data sets for which poor or incorrect results were obtained. There is still a lot to learn about the utility and limitations of the methods in these programs. - J E E, October, 31, 1989.