8.12 Spectral analysis, s(Spec)

The spectral analysis option for local and teleseismic events is selected in single trace mode. The spectral analysis is based on the Brune (1970) model and various assumptions about the geometrical spreading and anelastic attenuation. The steps in the analysis is:

$\bgroup\color{black}$\displaystyle d(f) = G(r,h) * D(f) * Moment*KK /(1+f**2/f0**2)* (4 * \pi * DE * V**3)) $\egroup$

where G(r,h) is geometrical spreading, r is epicentral distance, h is hypocentral depth, D(f) the diminution functiondue to anelastic attenuation, f is the frequency, DE the density, V the velocity at the source, f0 the corner frequency and KK a factor of 2.0*0.6 to correct for the free surface effect and radiation pattern.

$\bgroup\color{black}$\displaystyle D(f) = P(f) * exp (-\pi*f*trtime/(q0*f**qalpha)) where $\egroup$

trtime is the travel time from the origin time to the start of the spectral window and

is meant to account for near surface losses (Singh et al., 1982) with the constant kappa having a value of the order 0.02 sec. Anelastic attenuation Q is assumed to be frequency dependent following the relation $\bgroup\color{black}$Q = q0* f**qalpha$\egroup$ . For f less than 1 Hz, Q is, by default, assumed frequncy independent with the value q0 (before Dec 2013 it was frequecy dependent for f less than 1 Hz which could lead to a small overestimation of seismic moments for events larger than 4.5. This change is also affects other programs using Q-correcion like SPEC, AUTOMAG and AUTOSIG).

Before version 12.1, the transion of Q from high to low frequencies as a function of frequecy was made with the function Q=Q0*(1+f)**qalpha. However, this has turned out not to be ideal since Q at 1Hz then would be 2Q0 althought following the the function Q=Q0*(1+f/X)**qalpha at higher frequecies and becoming contant at lower frequencies. So now there is no smooth function and Q=q0 below 1 Hz and Q0 at 1 Hz. This parameter is set in SEISAN.DEF . For teleseismic events, only t* is used and Q must be set to zero (not used). The t* parameter is then entered using the kappa variable and is usually set to 1.0 (same value is used for P and S). So kappa is no used for teleseismic events. The geometrical spreading has been defined to be dependent on the wave type with several possibilities, all made equivalent to a distance called geo_distance (GD) such that geometrical spreading is expressed as 1/GD. There are several possibilities for GD:

P-waves:
GD is the hypocentral distance $\bgroup\color{black}$(HD) = sqrt (r*r +h*h)$\egroup$ so body wave spreading is assumed.

S-waves:
The geometrical spreading has been made dependent on distance and depth. At short distances, the geometrical spreading is assumed to be body wave spreading. For distances beyond the Herrmann-Kijko distance (default of 100 km) and a shallow focus, the following relation is used:

$\bgroup\color{black}$\displaystyle G(r,h) = 1/sqrt(100*r) =1/GD for r > 100 km $\egroup$

which is commonly used (Herrmann and Kijko, 1983; Herrmann, 1985). This relation assumes surface wave dispersion for epicentral distances larger than 100 km. In SEISAN 100 km is the default, however it can also be set to any other value by the parameter HERKIJ_DISTANCE (see later).

The above relation breaks down if the depth is large or comparable to the epicentral distance and in that case body wave spreading is again assumed. In order to get a smooth transition from surface wave to body wave spreading, it is assumed that the relation changes nearly linearly from surface wave spreading to body wave spreading between the depths GEO_DEPTH1 to GEO_DEPTH2. For depth less than GEO_DEPTH1(default 50 km), Herrmann-Kijko spreading is assumed, for depths larger than GEO_DEPTH2 (default 100 km), body wave spreading is assumed with the transition in between. In each case the geometrical spreading term is given as the equivalent GD, which is also recorded in the database. These 3 parameters can be used to change geometrical spreading. If e.g. HERKIJ_DISTANCE is 10 000 km, body wave spreading is always used. For more info, see (Havskov and Ottemöller, 2010).

$\bgroup\color{black}$\displaystyle G(r)=1/GD where GD =(27 + dist)/0.0048 $\egroup$

where dist is epicentral distance in degrees. This approximation is only valid for depth is less then 100 km and dist larger than 30 degrees.

The are 2 possibilities, use the values in MULPLT.DEF or use a spectral model in SEISAN.DEF.

Use the values given in MULPLT.DEF: Since only one value can be given for attenuation (Q and kappa), the same attenuation is used for P and S. Likewise, if the hypocenter changes depth, it will be the same fixed values for all attenuation and velocity parameters. For a particular fixed situation e.g. where the user always analyze shallow events and always is using the same type of spectrum (usually S), this is adequate. However in a situation with large depth variations, it is desirable to change all parameters as a function of depth. It might also be useful to change attenuation as a function of phase type although there often is little difference between P and S attenuation.

Use values in SEISAN.DEF: The spectral model is defined in SEISAN.DEF. In order to activate the model, set parameter SPECTRAL MODEL in MULPLT.DEF to 1.0, then the spectral model is used instead of the parameters in MULPLT.DEF. The model gives velocities (P and S), attenuation (P and S) at different depths as will as kappa for P and kappa for S and the behavior of Q below 1Hz (P and S).

The Q below 1 Hz parameter can now be set to any frequecy. If set to 1.0, the Q function is as above when using MULPLT.DEF parameters but if the frequecy is set to another valus, it is then the frequency f1 below which Q becomes constant with the value Q(f1). If the Q below 1Hz parameter is set to zero, the Q functions will not become constant at any frequecy..

Remember that the attenuation given at a particular depth is NOT the attenuation at that depth but the average between source and receiver. Values used will be interpolated between values in model and the values used will be show in the plot and stored in the S-file. If the spec model is used, it will be indicated on the plot (lower right hand corner). if a spec model is not there, MULPLT.DEF values will be used and also indicated on the plot.

From the spectral parameters, source radius and stress drop can be calculated as follows:

$\bgroup\color{black}$\displaystyle Stress drop = 0.44 * Moment /(source radius)**3 $\egroup$

The spectral analysis is used in two ways. The first and most common is to make the attenuation and instrument corrected displacement spectrum and determine the flat spectral level OM0, and corner frequency f0 from which the seismic moment, source radius and stress drop can be calculated. The second option is to display the instrument corrected spectrum (displacement, velocity or acceleration) and model the spectrum for corner frequency and attenuation parameters. In this case no correction for attenuation should be made.

Spectral analysis to determine moment, source radius and stress drop:
Select the spectral option, s(Spec). Before the spectrum comes up, you will get a question of the type of spectrum wanted and a few other options.

Displacment: d Displacement spectrum
Velocity: v Velocity spectrum
Acceleration: a Aceeleration spectrum
Raw spectrum: r Spectrum without instrument correction
Change velocity: c Use other velocity for calcualation
Change density: e Use other density for calculation
Change moment: m Change moment (see modeling)
Noise Pow. spec: n Noise powe spectrum (see later)
Lin axis: l Change to linear axis
New spec f-lim: f Change frequecy limits for spectrum
Change spectrum: t Change to P or S spectrum if auto does not work
Autofit spectrum: s Make automatic spec fit

The most used possibilities are displacement (d), velocity (v) or acceleration (a). For determination of Moment etc, the displacement spectrum or the autofit option (see below) MUST be selected. Unless raw spectrum is selected, the spectrum will be instrument corrected. If no response file is available in CAL, a message will be displayed on the screen and the raw spectrum calculated. At this stage it is also possible to change the velocity from the MULPLT.DEF value or the moment to be used for spectral fitting as given in the S-file (see spectral fitting below). The velocity change is only affective for the current spectrum since the user might make both P and S spectra and it is up to the user to make sure that it is appropriate for P or S, however, the value is saved in the S-file. The density change is kept for the whole MULPLT session since there is no change with velocity and its value is also saved in the S-file. The spectrum shown will normally show both the spectrum from the selected time window as well as a noise spectrum from an identical length time window at the start of the trace. IF NO NOISE SPECTRUM is desired, select spectrum with capital S instead of s.

If auto spectrum is chosen (s), an attempt will be made to fit the Brune spectrum to the observed spectrum. The program will first try to find the frequency range of acceptable S/N and then find the best fit by grid search and the fit will be shown on the spectrum as well as the frequency range used. The automatic frequency range determination can fail, particularly for small events so a check should be made if the spectrum fit looks 'reasonable' and that the frequency range is ok. Choosing a slightly different window might fix the problem. The initial frequency range is limited by the parameter SPECTRAL F-BAND set in MULPLT.DEF so changing this parameter might also help. Note that the automatic frequency range determination will fail if the noise window before the P is not as long as the analysis window chosen. Automtic spectra can also be done for the whole event with command automag in EEV, AUTOMAG and AUTO from the prompt line.

com_spec.out: The complex spectrum with some additional information needed for surface wave analysis, must be displacement spectrum.
amp_spec.out : The real spectrum given as frequencies and amplitudes. The files are only generated if parameter SPECTRAL_OUTPUT is set in MULPLT.DEF. Setting this parameter will also generate an ASCII waveform file with the input signal used.

Power spectra: The above spectra can also be displayed as power spectra if capital letters are used. Using e.g. 'V' instead of 'v' will show the power velocity spectrum.

When the spectrum comes up (see example in Figure 8.9 , the axis units are log amplitude in nanometers-sec (displacement) versus log frequency (Hz). The cursor can be used to select the level, corner frequency and slope by defining the spectrum with a 3 point selection. This 3-point selection is finished with f, q or r with the same meaning as in picking mode. The spectral values are displayed on the screen once q, f or r is pressed. The abbreviations are

NOTE: The veleocity is the velocity at the source, so if the dataset contains both shallow and deep earthquakes, a single velocity will be an approximation if the spec model is not used. There are two solutions to this problem if the values in MULPLT.DEF are used. A: Use different MULPLT.DEF for deep and shallow events. Since the the MULPLT.DEF used is taken from working directory if available, you can have different directories when working with different depth earthquakes. B: Use the general MULPLT.DEF in DAT and change the velocity after making the spectum in the S-file. At the next update, the moment etc will be recalculated with the new velocity. On top of the general parameter is indicated which kind of spectrum is assumed, P or S. In order for the program to automatically determine which kind of spectrum to assume, there must be a P or S reading displayed on the screen near the time window analyzed. The reading must be within 10 sec of the start of the window. If both a P and S-reading is within 10 secs, the nearest phase is chosen. If it cannot be determined which kind of phase is analyzed, the user will get a question to select type of phase (can also be changed later when spectral choices come up) The determination of which phase influences the further calculation of geometrical spreading and moment (uses P or S-velocity).

If f is selected, the spectral values together with calculated moment etc are stored in the S-file at the next key press (see parameters below). Spectral values in S-files accumulate, since no old values are deleted !!!. This is because the spectrum might be made under different conditions (start time, time window etc). The input parameters for the spectral analysis is given in file MULPLT.DEF, which can be in either DAT or the working directory, see below. Additional parameters for geometrical spreading are given in SEISAN.DEF in DAT.

$\bgroup\color{black}$\displaystyle Moment = 4 * pi * DE * V**3 * 10**OM /( G(r,h) * KK) $\egroup$

where V is the seismic wave velocity at the source (P or S if P or S-spectrum respectively) and OM the spectral flat level on the attenuation corrected displacement spectrum.

$\bgroup\color{black}$\displaystyle Moment magnitude = 2/3 * log10(moment) - 6.06 $\egroup$

$\bgroup\color{black}$\displaystyle Moment magnitude = 2/3 * log10(moment) -10.73 $\egroup$

The moment is calculated in Nm, the source radius in km and the stress drop in bars. All results are written to the S-file. Below is an example if Nordic format is used:

Note that no special line has been created in the Nordic format. Comment lines are used with SPEC at the start of the line followed by station and component. Only the first 4 characters of the 5 character station name is used. The station and component names are given at the start of the line. In case of a 5 character station name, the station name is shifted one character to the left. The information is:

Note: Before version 10.1, VS was written instead of Q1. However, VS was not used for anything. Note: The component codes have not been adjusted for SEED so the location code is not included.

Note: In earlier versions (before version 7.0), the field for kappa was used for the travel time to start of window. This can be calculated from origin time and the start time of the window.

In Nordic2 format, a special spectral line (type S) has been made and now there is also room for network and location. Example below:

Now the spectral information is on one line and there is an explanatory line above with the same abbreviations as in Nordic format except VX which now indicates P or S velocity. SEISAN can read both new and old format independent of whether Nordic or Nordic2 format is used but will only write S-lines if Nordic2 format is used. This should ensure backwards compability.

When doing an UPDATE of the database or just a location with HYP, all distance dependent spectral values are recalculated and average values written into the output file. Mw will be calculated from the average value and written in the header line. However, the original distance dependent Q and kappa correction is not changed, since this correction was used to modify the spectrum used for reading parameters. Normally a small distance change has insignificant influence on the spectral level or the corner frequency so the Q-correction should be no problem. Spectra of the same type (P, S or ?) and from the same channel are overwritten. Only in case of UPDATE are the values written back into the database.

Display of spectral parameters: Program MAG can read and plot relations between spectral and source parameters. Program REPORT can read spectral parameters and combine in a table.

Potential problem with Q-correction: If the origin time in header is wrong, the Q-correction can be very wrong.

There must be a phase line in the S-file with component and distance corresponding to the spectra made in order for the spectral values to be calculated.

Once the spectrum has been shown (displacement, velocity or acceleration), a theoretical spectrum can be calculated and superimposed on the observed spectrum in order to forward model either source parameters or attenuation.

Another option (s: Plot stored fitted spectrum) if spectral analysis already has been done, is to plot the corresonding theoretical spectrum on top of the calculated spectrum. This option has the purpose to be able to check previously made spectra whether manual or automatic. The SPEC lines have all the information needed to make the theortical spectrum.

Entering constants and modeling: The modeling can only take place when the spectrum is seen on the screen.

Press S and a question will appear to enter the constants f0, k, Q0 and qa which are as defined above except qa is Qalpha. Once these parameters have been entered (terminate with return), the theoretical spectrum (displacement, velocity or acceleration depending on what is used for the spectrum) is calculated and superimposed on the observed spectrum. The parameters used or calculated are displayed. S or s can now be pressed again and a new theoretical spectrum calculated and plotted. To get out of the spectral fitting loop, type r or q as usual.

Which constants and parameters are used: The moment is taken from the last moment calculation made for thsi station, if any. If none, the moment is taken from the S-file if an average moment has been calculated (see UPDATE command). If no moment is available, it can also be entered the first time the spectrum is shown (option M). If no moment is available, no modelling can be done and a message is given. The distance and depth is likewise taken from the S-file. If no distance is available, no modelling can be done and a message is given. If all 4 parameters f0, k, Q, qa are entered, stress drop is calculated with the relation given above. If the corner frequency is given as zero, the user will be asked to enter the stress drop and the corner frequency is calculated from the stress drop. If Q is zero, no Q-correction is made. IMPORTANT: The Q and qa used here are distinct from the Q0 and Qalpha used for making the amplitude spectrum and both should not be used when modeling since this would imply a Q-correction two times. The best way is to use Q0=0 and kappa=0, so that Q is only corrected for when modeling. The distance used is everywhere is GEO_DISTANCE.

The 3 types of spectra (displacement, velocity and acceleration) can optionally be made as power spectra. Instead of selecting the type of spectrum by pressing d, v or a, just press the same characters in upper case and the power spectrum will be shown.

In seismic noise studies, the seismic background noise is often displayed as acceleration power spectral density in dB relative to ((1m/s**2)**2)/Hz. Instead of selecting d, v or a, press n instead. The plot shows the Peterson (1993) new global high and low noise models superimposed on the observed spectrum (Figure 8.8). When doing noise spectra, no attenuation correction is done. The normalization of the spectrum is as follows

$\bgroup\color{black}$\displaystyle P = \vert F^{DFT} \vert^{2} \times \frac{\Delta t^{2}}{T} \times 2 $\egroup$

where P is the Peterson Power spectrum, $\bgroup\color{black}$F^{DFT}$\egroup$ is the discrete Fourier transform, $\bgroup\color{black}$\Delta t$\egroup$ is the sample interval and T is the length of the time window. The factor 2 comes from the fact that only the positive frequencies are used so only half the energy is accounted for. The total power is proportional to the length of the time window since the noise is considered stationary, so by normalizing by T, the length of the time window should not influence the results. This noise option is a handy method of checking the noise characteristics of a given seismic station and compare it to global standards. This kind of analysis can also be done with the SPEC program (section 32). For more information, see instrument.pdf in INF.

**Figure 8.8:** Example of a noise spectrum.
$\begin{figure} \centerline{\includegraphics[width=0.9\linewidth]{fig/fig23}} \end{figure}$

Problems: There is currently no check if a displacement seismogram has been calculated when calculating the spectral parameters. If spectral analysis is done outside EEV (output in MULPLT.OUT) or with EEV when there is no origin time and/or epicentral distance, the output results are wrong for moment etc. Before calculating moment etc, the S-file MUST HAVE BEEN UPDATED SINCE BOTH THE DISTANCE AND ORIGIN TIMES ARE USED. If the spectra get very high amplitude levels when correcting for instrument, this might be caused by correcting for Q. With a Q of 100 and a distance of 10 000 km, this gives a very large correction. The Q-correction can be disabled in the MULPLT.DEF file. If picks are made, but no readings appear in the S-file or readings appear with wrong component, the waveform file component might not have been defined in subroutine componen.for. If poles and zeros are used to remove the response, rotation cannot be used at the same time.

**Figure 8.9:** Spectral analysis
On top the original trace is seen and on the bottom the displacement spectrum (log -log, unit nm-sec and Hz). The level and slope has been indicated interactively. Note the noise spectrum at the bottom of the figure.
$\begin{figure} \centerline{\includegraphics[width=0.9\linewidth]{fig2/fig4a}} \end{figure}$

**Figure 8.10:** Three component analysis
On top the Z-channel is shown together with the window used for the 3 channels Z, N and E shown below. The signals below has been filtered between 1 and 5 Hz and the resulting azimuth of arrival is 160 degrees and a correlation coefficient of 0.2. The apparent velocity is 9.8 km/sec.
$\begin{figure} \centerline{\includegraphics[width=0.9\linewidth]{fig2/fig4b}} \end{figure}$

**Figure 8.11:** Plotting response curves
The figure shows the amplitude and phase response for station SUE, component S Z. The response is the one which will be used in analysis irrespective of whether it is taken from the file header or the CAL directory.
$\begin{figure} \centerline{\includegraphics[width=0.9\linewidth]{fig2/fig4c}} \end{figure}$

Vel:	Velocity used (km/sec) (Vp or Vs)
Dens:	Density (g/(cm**3)
Dist:	Hypocentral distance (km)
Depth	Hypocentral depth (km)
q0:	q0 for spectral amplitude correction
qalpha:	qalpha for spectral amplitude correction
k:	kappa
q<1Hz:	frequency where Q(f) changes to constant Q
spec model:	Spec model is used

MO:	log of moment, unit Newton*m
ST:	Stress drop in bars
OM:	log spectral level (nm-sec) not distance corrected
F0:	Corner frequency (Hz)
R :	Source radius (km)
AL:	Decay of log spectrum
WI:	Spectral window used (secs)
MW:	Moment magnitude
T :	Start time of window for spectrum in hr, min, sec
K:	Kappa
GD:	GEO_DISTANCE in km
VP or VS:	Velocity in km/sec at source for P and S-spectra respectively. The P or S in this line indicated if the spectrum is a P-spectrum or an S-spectrum. It MUST be P or S to be used for magnitude determination. A `?' is put in if MULPLT does not know which kind of spectrum is (no P or S reading near start of spectral window). This can be changed by editing the S-file afterwards.
DE:	Density in g/cm**3
Q0:	q0 in relation Q = q0 * f ** qalpha
QA:	qalpha
Q1:	frequency where Q(f) changes to constant Q

Moment:	Moment used
Geo dist:	Geo distance used
Stress drop:	Stress drop in bars
f0:	Corner frequency
k:	Constant used in diminution function
q:	q0 used in spectral fitting
qa:	qalpha used in spectral fitting