Calculating Auxiliary Fault Plane Solutions given the main fault (Fortran)

We have calculated the auxiliary fault plane solution using the input of main fault solution in perl. Here, we do the same in Fortran. And we will also plot to solution to visualize the results.

Fortran Code to get the auxiliary fault plane solutions:

program auxiliary_fault_plane
! Program to calculate the strike, dip and rake of the auxiliary fault plane solutions
! given the strike, dip and rake of the fault plane.
! Authors: Utpal Kumar, Li Zhao
        implicit none
        real (kind=8):: ph1,del1,lb1,ph2,del2,lb2,sinlb2,coslb2,sinph1_ph2,cosph1_ph2,ph1_ph2,phI,delI,lbI
    real, parameter :: pi = 2*asin(1.0)
    integer :: n,lbsgn

        print *, 'Enter strike, dip and rake respectively of the fault plane (in degrees) (e.g., 40,50,60):'
        read *, ph1,del1,lb1
        !ph1=40
        !del1=80
        !lb1=200
        phI=ph1
    delI=del1
    lbI=lb1
        !Converting the strike, dip and rake in radians
        ph1=ph1*pi/180  !strike
        del1=del1*pi/180        !dip    
        lb1=lb1*pi/180          !rake   
        
        !Adaptation to include the negative values of rake
        lbsgn=1
        if (lb1 < 0) then
                lbsgn=-1
        end if
        
        
        ! Calculation of Strike, sip and rake of the auxiliary fault plane
        del2=acos(sin(lb1)*sin(del1))   !dip of auxialiary fault plane
        sinlb2=cos(del1)/sin(del2)
        coslb2=-(sin(del1)*cos(lb1)/sin(del2))
        lb2=acos(coslb2)        !rake of auxilairy fault plane
        sinph1_ph2=cos(lb1)/sin(del2)
        cosph1_ph2=-1/(tan(del1)*tan(del2))
        ph1_ph2=acos(cosph1_ph2)

        !Checking for the quadrant of the strike angle
        if (sinph1_ph2 >= 0 .and. cosph1_ph2 >=0) then
        ph1_ph2=ph1_ph2
        else if (sinph1_ph2 > 0 .and. cosph1_ph2 < 0) then
        ph1_ph2=ph1_ph2
        else if (sinph1_ph2 < 0 .and. cosph1_ph2 < 0) then
        ph1_ph2=-ph1_ph2
        else if (sinph1_ph2 < 0 .and. cosph1_ph2 > 0) then
        ph1_ph2=-ph1_ph2
        end if
        
        ph2=ph1-ph1_ph2         !strike of auxialiary fault plane
        
        ! For dip > 90 degrees and less than 180 degrees
        if (del2 > pi/2 .and. del2 < pi) then
        ph2= pi + ph2
        del2= pi - del2
        lb2= 2*pi - lb2
        end if
        
        if (lbsgn < 0) then
                lb2 = -(2*pi - lb2)
        end if
        
        !Adaptation to give the strike value in the range of 0 to 360 degrees
        if (ph2 > 2*pi) then
                ph2 = ph2 - 2*pi
        end if
        
        101 format("The strike of the auxiliary plane is: ",f9.4, " degrees")
        102 format("The dip of the auxiliary plane is: ",f9.4, " degrees")
        103 format("The rake of the auxiliary plane is: ",f9.4, " degrees")
        print 101, ph2*180/pi
        print 102,  del2*180/pi
        print 103, lb2*180/pi
    open(unit=10,file='plt.dat')
    10 format("25  25   0 ",  f8.2, f8.2, f8.2, f8.2, f8.2, f8.2," MainFault: ",f7.2,"/",f7.2,"/",f7.2, &    !continuation of line
        " AuxFault: ",f7.2,"/",f7.2,"/",f7.2)
        write(10,10) phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi, phI,delI,lbI,ph2*180/pi, del2*180/pi, lb2*180/pi !writing in the file
end program auxiliary_fault_plane

Then we save this in the file called “auxiliary_fault_plane.f”. This script basically aim for calculation of the fault plane solution. For visualizing the results, we can execute the program and plot the results using the psmeca command of GMT.

The bash script for compiling and executing the above program and then plotting it is following:

#!/bin/bash
gfortran -ffree-form auxiliary_fault_plane.f -o auxiliary_fault_plane
./auxiliary_fault_plane
psbasemap -Bwsne -R0/50/0/50 -Jm1.0 -Xc -Yc -P -K >output.ps
echo "24.6 30 20 0 2 MC Focal Mechanism Plot"| pstext -Jm -R -K -O -P -N >>output.ps
psmeca plt.dat -Jm -R -Sa0.40/14/-0.2i -Gbrown -Fa0.2c/cd -Egray -P -O -N -L -V >>output.ps

gs output.ps
ps2pdf output.ps

The P axis is plotted using circle and T axis using the diamond symbol (-Fa0.2c/cd).

-Utpal Kumar (IES, Academia Sinica)

Using Generic Mapping Tool (GMT) for Oblique Projection Problem

In previous post (Using Generic Mapping Tool (GMT) Basic) we make a map with the study area is located at Taiwan, where the subduction trench are either having South-North and or West-East  orientations, respectively. For these subduction systems, we could make a cross section immediately perpendicular to the trench orientation, i.e. West-East (for the South-North oriented trench) and South-North (for the West-East orientated trench). From this cross section, we could see the Wadati-Benioff zone clearly from the earthquakes distribution (although clear or not will be depend on your data quality).

However, this not the case for the study area at Sumatra (Using Generic Mapping Tool (GMT) Basic II), where the subduction trench is Northeast-Southwest orientation or we can call it as an oblique subduction. If we make a cross section as same as the one that we described earlier, we wouldn’t able to see the Wadati-Benioff zone. Now this is the problem, how we could rotate the map so that it will have a parallel orientation to the trench?

For this problem, we could solve it easily by using the GMT. The GMT provide various map projection that enable us to plot a map in almost what ever orientation that we wanted to. Keep in mind that, in this post we work on a Linux Operating System (OS) with GMT version 4.5 (some features, i.e., to make a file executable, post script viewing (.ps) and others, you may need to change it or install some program in order to work and viewing the output).

For this post, the files that you need are listed as shown below (freely accessible).

1. Topography grid file GEBCO 30-arc second (download here).
2. Earthquake event (customized for this post only) at northern Sumatra region from the USGS earthquake catalog (download here).
3. Sumatra trench file from the USGS (download here).
4. Color file for the map (you can create by yourself, or download here).
5. Color file for the earthquake event (same as number 4, or download here).

After you downloaded all the files, make sure that these files are located in the same folder as your scrip. Now, you may try to plot the Map by using this script as follow.

#!/bin/csh
#environment setting
gmtset PAPER_MEDIA A3
gmtset LABEL_FONT_SIZE = 12p
gmtset BASEMAP_TYPE plain
gmtset ANNOT_FONT_SIZE_PRIMARY = 12p
gmtset ANNOT_FONT_SIZE_SECONDARY = 12p
gmtset HEADER_FONT_SIZE = 14p

####topographic file####
set inputgrd=sumatra1.grd

####range and projection for regional map####
set range1=95/-2/97/8r
set oblq1=98.5/4/55/4i

####range and projection for study area map####
set range2=96/-1/95.95/5r
set oblq2=100/10/10/54/5.8i

####cpt, gradient and trench####
set mapcpt=sumatra.cpt
set inputgradient=sum1_intens.grd
set trench=sumatra_trench.clip
set eventcpt=eventdepth.cpt

####erthquake event from usgs cataloque custom area and time#####
#magnitude 4
awk '{if ($5>=4 && $5<5) print $3, $2, $4}' usgs_cat2005.dat > m4_2005.xyz 
#magnitude 5
awk '{if ($5>=5 && $5<6) print $3, $2, $4}' usgs_cat2005.dat > m5_2005.xyz
#magnitude 6
awk '{if ($5>=6 && $5<6) print $3, $2, $4}' usgs_cat2005.dat > m6_2005.xyz
#magnitude 7
awk '{if ($5>=7 && $5<8) print $3, $2, $4}' usgs_cat2005.dat > m7_2005.xyz
#magnitude 8
awk '{if ($5>=8 && $5<9) print $3, $2, $4}' usgs_cat2005.dat > m8_2005.xyz

####output####
set output=oblique_subduction.ps

####Map Script for regional map####

#grd image#
grdgradient $inputgrd -A0/270 -G$inputgradient -Ne0.6 -V
grdimage $inputgrd -R$range1 -JOa$oblq1 -Ba1f0.5WNSE -C$mapcpt -I$inputgradient -X11.7i -Y7.5i -K -V > $output

#coast line#
pscoast -R -J -B -Dh -W1.5p,black -O -K -V >> $output

#custom contour line#
grdcontour sumatra1.grd -R -J -C500 -L-7000/-4500 -K -O -V >> $output
grdcontour sumatra1.grd -R -J -C100 -L-4500/-3500 -K -O -V >> $output

#sumatra_trench#
psxy $trench -R -B -J -Sf0.5i/0.1ilt -Gred -W0.03i,red -K -O -V >> $output

#input event custom sized#
#magnitude 4
psxy m4_2005.xyz -R -J -Sc0.04i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 5
psxy m5_2005.xyz -R -J -Sc0.07i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 6
psxy m6_2005.xyz -R -J -Sc0.10i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 7
psxy m7_2005.xyz -R -J -Sc0.15i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 8
psxy m7_2005.xyz -R -J -Sc0.20i -W0.001p -C$eventcpt -K -O -V >> $output

#make a rectangle indicate for the study area#
psxy -R -J -W0.02i/red -K -O -V << EOF >> $output
95.95 5
93.20 2.98
96.1 -1
#> next (this gonna disconnect the points)#
98.8 1.01
95.95 5
EOF

#cross section line#
psxy -R -J -W0.02i,black -K -O -V << END >> $output
94.6 1
97.4 3
END

####Map Script for study area map####
grdimage $inputgrd -R$range2 -JOc$oblq2 -Ba1f0.5WNSE -C$mapcpt -I$inputgradient -X-9i -Y-4.3i -K -O -V >> $output

#coast line#
pscoast -R -J -Dh -Ba1f0.5WNSE -W1.5p,black -K -O -V >> $output

#custom contour line#
grdcontour sumatra1.grd -R -J -C500 -L-5000/-4500 -K -O -V >> $output
grdcontour sumatra1.grd -R -J -C100 -L-4000/-1000 -K -O -V >> $output
grdcontour sumatra1.grd -R -J -C250 -L-1000/-250 -K -O -V >> $output

#input sumatra_trench
psxy $trench -R -B -J -Sf0.5i/0.1ilt -Gred -W0.03i,red -K -O -V >> $output

#input event custom sized#
#magnitude 4
psxy m4_2005.xyz -R -J -Sc0.04i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 5
psxy m5_2005.xyz -R -J -Sc0.07i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 6
psxy m6_2005.xyz -R -J -Sc0.10i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 7
psxy m7_2005.xyz -R -J -Sc0.15i -W0.001p -C$eventcpt -K -O -V >> $output
#magnitude 8
psxy m8_2005.xyz -R -J -Sc0.20i -W0.001p -C$eventcpt -K -O -V >> $output

#cross section line#
psxy -R -J -W0.02i,black -O -V << END >> $output
94.6 1
97.4 3
END

gv $output

After you execute the program, you may get a map as shown below.

Here are some explanations to the script:

1. Oblique projection can be made by using range as described in set range1 and set range2, which are the range of the regional and study area maps respectively. If you see carefully, we append a small “r” later at the end of the range. This implies that we will use lower left corner and upper right corner longitude and latitude information for the map range.
2. We use two types of oblique projection (you may type in the terminal man grdimage to get a full description):
   1. “-JOa$oblq1” means that we will use a cylindrical projection with an Oblique Mercator – point and azimuth.
   2.  “-JOc$oblq2” same as 1, but with an Oblique Mercator – point and pole.
3. We make the illuminations of the topography grid file by using grdgradient.
4. We use a simple shell command such as awk, i.e. to get the preferred column being printed to the file that we want, such as the usgs_cat2005.dat being printed to a several files with different magnitude (column 5 in the .dat file indicate for the earthquake magnitude).
5. We set up a custom earthquake size. In default, the psxy command will search the 4th column as the earthquake size and will plot the size using default size. However, we can customized the size that we want to see the earthquake with larger magnitude has a “logarithmic” size from the small one, so it will give us an “appealing” attention to the big magnitude event as the main shock.

Keep in mind that, this post only show how to rotate the map, not the cross section plot. However, we already make a line for the cross section and we will cover that for the next post.

•                                  Haekal A. Haridhi, TIGP Program, IES, Academia Sinica

Analytical signal and Hilbert Transform

For a real time series x(t), its analytic signal x(t) is defined as

x(t) = x(t) – iH[x(t)]

Let us consider an example of a monochromatic signal 𝑥(𝑡) = 5 sin(10𝑡 + 3).

Figure 1

Now, let us consider a more complex function

x(t) = 1*sin(2𝜋10𝑡 + 0.3) + 2*sin(2𝜋20𝑡 + 0.2) + 3*sin(2𝜋30𝑡 + 0.4).

Figure 2

We can clearly observe that the Hilbert transform estimates the instantaneous frequency of a signal for monocomponent signals only.

Matlab Codes:

clear; close all; clc

fs = 1e4;

t = 0:1/fs:1;

f1=10;

f2=20;

f3=30;

x=1*sin(2*pi*f1*t+0.3)+2*sin(2*pi*f2*t+0.2)+3*sin(2*pi*f3*t+0.4);

%x=5*sin(2*pi*10*t+3);

y = hilbert(x);

figure(1)

plot(t,real(y),t,imag(y))

%xlim([0.01 0.03])

legend('real','imaginary')

title('Hilbert Function')

figure(2)

subplot(3,1,1)

env=abs(y);

plot(t,x)

xlabel('Time')

title('Envelope')

hold on

plot(t,env)

legend('original','envelope')

subplot(312)

instph=fs/(2*pi)*unwrap(angle(y));

plot(t,instph)

xlabel('Time')

ylabel('Phase (in rad)')

grid on

title('Instantaneous Phase')

subplot(313)

instfreq = fs/(2*pi)*diff(unwrap(angle(y)));

plot(t(2:end),instfreq)

xlabel('Time')

ylabel('Hz')

grid on

title('Instantaneous Frequency')

Application of Hilbert Transform:

Figure 3

Aliasing and Nyquist Condition

Let us consider a simple sinusoidal signal

x(t) = 5sin(10t + 3)

This is a signal with an amplitude of 5, the angular frequency of 10 and phase of 3. The Nyquist condition demands that the angular frequency, mod(𝜔) ≤ pi/delta t.

For the above figure 1, we took 𝜔 to be 10 and pi/delta t is 99.5, so the Nyquist condition satisfies and we do not see the aliasing. Here the sampling rate, Δ𝑡 is 0.0316.

Let us see what happens when we take the pi/delta t to be less than 𝜔.

In this case, we took the sampling rate, Δ𝑡 to be 0.3307. Here the pi/delta t is 9.5000 which is less than 𝜔. Hence the Nyquist condition is not satisfied and there is a clear aliasing. The sinusoidal signal does not appear to be sinusoidal.

Figures

Singular Value Decomposition

The inversion of a matrix is a very important part of data analysis. The quality of inverse solution depends on how well we can invert the matrix.

d=Gm where “d” is the N dimensional column vector consisting the data, “m” is the M dimensional column vector consisting the model parameters which we seek to invert for, “G” is the NxM dimensional kernel matrix to map the model vectors into the data matrix.

=>m = Ad, where A is the inverse of the matrix G.

There are several matrix decomposition techniques. These techniques are not only useful for the inverse problems, but also for applications like principal decomposition that are widely used in seismic attribute studies.

When a matrix is pre-multiplied to a vector, the resultant vector can be regarded as a linear transformed version of the original vector. Therefore the multiplication of matrices is a linear transformation. Any matrix can generally be decomposed into transformations of some other matrices.

The decomposition of a general rectangular matrix requires the use of singular value decomposition or SVD (Lanczos, 1961). The SVD decomposes any rectangular matrix A of m rows and n columns into a multiplication of three matrices of useful properties:

A = U L V’

Here, L is the m x n rectangular diagonal matrix containing p singular values (or principal values) of the matrix A.

U and V are the two square unitary matrices. The number of non-trivial singular values, p is called the trace or rank of the matrix.

We try to solve for the an example of

A = (10  2; -10 2)

Example 1

Understanding Seismograms

Seismograms are basic information about earthquakes, chemical and nuclear explosions, mining induced earthquakes, rock bursts and other events generating seismic waves.

Seismograms reflect the combined influence of the seismic source, the propagation paths, the frequency response of the recording instrument and the ambient noise at the recording site.

u(t) = s(t)*g(t)*i(t)*n(t).         “*” means convolution.

where, u is the seismic record, s is the source effect, g is the propagation effect, i is the instrument effect and n is the noise.

So, by understanding seismogram we can understand the seismic source, earth structure, and the noise in the medium.

This is our short attempt to understand the seismograms.

Understanding Seismograms 1

For understanding the seismogram, there are two main things we should notice- record duration and the dispersion. Due to different nature of the propagation velocity of the seismic waves and the different propagation paths taken by them to the station, travel time differences between the main group usually grows with distance. Since the body wave groups do not disperses, so their individual duration remains more or less constant, only the time difference between them changes with distance.

Figure 1: Event-Station Distribution map for the event in Mindanao, Phillipines Island (2005-02-05, Mw-7). The source is 540 km deep.

Figure 2

Figure 3

Figure 4

The epicentral distance for QIZ, BJT and DGAR is ~ 19, 35, 52 degrees respectively. For the station BJT (Figure 3), the epicentral distance is 35.2 degrees and for the station DGAR (Figure 4), the epicentral distance is 52.4 degrees. The P-S arrival times for the station BJT is nearly 300 sec and for the station DGAR, the P-S arrival time is nearly 400 sec. For the station QIZ, it is only around 200 sec. So, with the increase in epicentral distance, the difference in individual phase arrival increases. Also for the station DGAR, the P arrival is later than the BJT.

The time difference between main body wave onsets is:

Distance (degrees) less than	Time-difference in body wave onsets (P-S) in sec
10	180
60 (Fig 1,2,3)	960
100	1800
180	2700

Figure 5 : Travel-Time curve for depth 540 km

For the deep events, we do not see the first arrival as P wave (Figure 6) for short distances (D

Figure 6

Figure 6, is the seismogram at 2 degrees epicentral distance from the deep source. Here, we can notice that the first arrival is not P but the core reflected phase PcP. It is also cevident from the travel time curve(Figure 5).

Figure 7: Event-Station Distribution Map of Near Coast of Peru (2001-06-23, Mw-8.3). The source is at 2 km depth.

When the source is near the surface then we observe the surface waves as well in the seismogram.The surface waves from earthquakes at intermediate (> 70 km) or great depth (> 300 km) may have amplitudes smaller than those of body waves or may not even be detected on seismic records.

In contrast to the body waves, the velocity of surface waves is frequency dependent and hence it is dispersed. The duration of Love and Rayleigh waves increases with distance. As for station NNA (Figure 8), the duration of the first group of Rayleigh wave is about 230 seconds and for station HKT (Figure 9), the duration is around 600 seconds.

We can also notice that for shallower earthquake event, P-wave is the first arrival for stations less than 10 degrees of the source (Figure 10).

Figure 8

In Figure 8, we can notice that for surface waves, the largers periods arrive first than the lowers periods. This is the general case for normal layering. The larger periods sample the higher depths and since the velocity for higher depths is larger so it arrives before the shorter period waves.

Figure 9

In figure 9, for the surface waves, we can observe the wave groups or the beating effect (analogous to acoustics). This is because of the interference of the two or more harmonic waves with closer frequencies or periods.

Figure 10: Travel time curve for depth 2.2km

For certain distance ranges, the travel time curves for different types of seismic phases are close to each other or can even overlap (Figure 5,11).

(a)

(b)

 Figure 11: S, SKS and ScS are very close to each other.

The amplitude of P-wave depends not only on the epicentral distance but on the site effect (underground geology and crustal heterogeneity), azimuth etc.

The most obvious indication on a seismogram that a large earthquake has a deep focus is the small amplitude of the surface waves with respect to the body-wave amplitudes and the P and S waveforms often have impulsive onsets.

A more precise determination of the depth h of a seismic source, however, requires either the availability of a seismic network with at least one station being very near to the source, e.g., at an epicentral distance D < h (because only in the near range the travel time t(D, h) of the direct P wave varies strongly with source depth h), or the identification of seismic depthphases on the seismic record.

At distant seismograph stations, the depth phases pP or sP follow the direct P wave by a time interval that changes only slowly with distance but rapidly with depth.

Inverse Theory

We follow the Geophysical Data Analysis: Discrete Inverse Theory by William Menke book.

We solve the chapters and try to make the things more simpler and faster to read and understand.

Chapter 1: Probability Theory

Click the link below to download the notes:

Probability Theory

Chapter 2: Solution of the Linear, Gaussian Inverse Problem : The length method

Linear, Gaussian Inverse Problem 1

Chapter 3: Solution of the Linear, Gaussian Inverse Problem : Generalized Inverses

Linear, Gaussian Inverse Problem 2

Stay tuned for the other chapters…

Moment Tensor Calculations

Programs:

1. moment_tensor_calc2.pl –> it calculates the moment tensor matrix for given M0, strike, dip and rake

2. diagonalizing_matrix.pl –> this program diagonalizes the given matrix. It can be used to diagonalize the moment tensor matrix as well.

3. iso_dev_mt.pl –> it gives the isotropic and deviatoric part for a given moment tensor matrix

4. maj_min_dc.pl –> this gives the major and minor double couple as output.

5. dc_clvd_mt.pl –> this gives the double couple and clvd part as output.

6. focal_mech_info.pl –> this gives the auxiliary fault plane as output for a given main fault plane.

7. plot_focal_mech.sh –> this bash script plots the focal mechanism solution for a given strike, dip and rake. It makes use of the focal_mech_info.pl program.

Example Focal Mechanism Plot for:

Main Fault

Strike = 40 degrees

Dip = 50 degrees

Rake = 60 degrees

Click here for the plot