Helmer Cluster show case 3: 24 micro earthquake released energy within 100 ms around a tunnel

The analytic solutions do not exist in the modeling situation where a tunnel dipping toward the east and ends somewhere within the model (See Fig. 1a). The total CPU time on my 18 core cluster is about 15 hours.

Fig1a: Tunnel Model with acquisition geometry

Fig1a: Tunnel Model with acquisition geometry

Fig. 1b provides a closer view of the tunnel and the spatial distribution of the 24 micro-earthquake sources. The “delay” on the figure is the original time delayed with respect to the first source. Source types are: ISO-explosion, DC-double couple, X-linear vector along X direction, and Z-linear vector along vertical direction. The frequencies ranges from 50 Hz to 300 Hz.

Fig1b: A close view of the tunnel and the spatial distribution of the micro-earthquakes

Fig1b: A close view of the tunnel and the spatial distribution of the micro-earthquakes

Fig. 1c and 1d shows the two components seismic traces acquired by the surface receivers and the borehole receivers, respectively. It is easy to observe that the second phase immediately following the first P-wave phase is not necessarily the shear-wave phase of the same earthquake. Instead, it can be a P-wave phase of a late earthquake.

Fig1c: Seismic gather from the surface receiver array. Fig 1d: Seismic gather from the borehole receiver array.

Fig1c: Seismic gather from the surface receiver array. Fig 1d: Seismic gather from the borehole receiver array.


Fig 1e and 1f shows the waveform from the middle channel from the surface receiver array (channel 76) and from the borehole receiver array ( channel 51). Can you identify all the phases in this noise-free environment? I actually lied when I say I used 24 sources. Can you find out how many sources I really used?
Fig1e: The waveform of the two components surface receiver. Fig 1f: The waveform of the two components borehole receiver.

Fig1e: The waveform of the two components surface receiver. Fig 1f: The waveform of the two components borehole receiver.

Here is a movie showing you what is really going on inside the earth, which in reality you never observe.

Still using old-fashion Hudson TK plot?

Hudson TK plot (fig. 1) was proposed by Hudson, Pearce, and Rogers, in 1989 to represent moment tensor patterns on a 2D uv-plane. T and k, parameters describing the moment tensor pattern can be calculated from uv. Recently, Tape and Tape (2012) proposed that displaying patterns on a sphere is more natural. I agree. Plotting moment tensor patterns on a lune surface on a sphere is more intuitive than on a TK-plot which actually uses uv as coordinate and the connection of TK to the eigenvalues of the moment tensor is less straightforward.

Fig 1. TK-plot

From now on, I will call the conventional Hudson method as HPR method and the lune surface method as RJ method, since the spherical representation was first used by Riedesel and Jordan (1989). I write a matlab script that is ready to generate a lune mesh on a sphere. See fig. 2.

Fig. 2: lune mesh on a sphere with typical moment tensor marked on the surface, i.e., Isotropic (ISO), LVD (linear vector), CLVD (compensated linear vector dipole).

Matlab Script:

function fig=lunemesh
RR=inv([sqrt(3) 0 -sqrt(3);-1 2 -1; sqrt(2) sqrt(2) sqrt(2)]/sqrt(6));

fig =figure;set(fig,'units','normalized','outerposition',[0 0 1 1]);
[Xout, Yout, Zout, ~]=sphere3d(ones(37),-pi,pi,-pi/2,pi/2,1,1,'off','spline',0);
tmp=RR*[Xout(:)';Yout(:)';Zout(:)'];
Xout(:)=tmp(1,:);Yout(:)=tmp(2,:);Zout(:)=tmp(3,:);
mesh(Xout,Yout,Zout,zeros(size(Xout)));
hold on
Xl=zeros(size(Xout));Yl=Xl;Zl=Xl;

Xl(:,16:22)=Xout(:,16:22);Yl(:,16:22)=Yout(:,16:22);Zl(:,16:22)=Zout(:,16:22);

mesh(Xl,Yl,Zl,ones(size(Xl)),'linewidth',1);
colormap ([.8 .8 .8;0 0 0]);

% add arrow
[hline_x,hhead_x]=arrow3d([0 0 0],[2 0 0],15,'cylinder',[0.1,0.15]);
set(hhead_x,'facecolor','red','edgecolor','black');set(hline_x,'facecolor','black');
[hline_y,hhead_y]=arrow3d([0 0 0],[0 2 0],15,'cylinder',[0.1,0.15]);
set(hhead_y,'facecolor','blue','edgecolor','black');set(hline_y,'facecolor','black');
[hline_z,hhead_z]=arrow3d([0 0 0],[0 0 2],15,'cylinder',[0.1,0.15]);
set(hhead_z,'facecolor','yellow','edgecolor','black');set(hline_z,'facecolor','black');

axis auto
axis equal
axis([-1 1.5 -1 1.5 -0.2 0.1])
axis off

caxis([0 1])
text(2.2, 0, 0,'lambda_1','fontsize',15)
text(0, 2.2, 0,'lambda_2','fontsize',15)
text(0.3, 0.0,1.8, 'lambda_3','fontsize',15)

campos([ 17.7381    7.3697  -17.9197]);
camup([ 0.4829    0.5507    0.6808])
camzoom(1.2)
% add big arc

ISO1=[1,1,1];ISO2=[-1 -1 -1];
ISO1=ISO1/norm(ISO1); ISO2=ISO2/norm(ISO2);

C1=[3,1,1];C2=[-1 -1 -3];
C1=C1/norm(C1);C2=C2/norm(C2);

LVD1=[1,0,0];LVD2=[0 0 -1];
LVD1=LVD1/norm(LVD1); LVD2=LVD2/norm(LVD2);

CLVD1=[2 -1 -1];CLVD2=[1,1,-2];
CLVD1=CLVD1/norm(CLVD1);CLVD2=CLVD2/norm(CLVD2);

Carc=vec2arc(C1,C2);
LVDarc=vec2arc(LVD1,LVD2);
CLVDarc=vec2arc(CLVD1,CLVD2);

plot3(Carc(:,1),Carc(:,2),Carc(:,3),'k.-','linewidth',1.5);
plot3(LVDarc(:,1),LVDarc(:,2),LVDarc(:,3),'r.-','linewidth',1.5);
plot3(CLVDarc(:,1),CLVDarc(:,2),CLVDarc(:,3),'c.-','linewidth',1.5);

plot3(ISO1(1),ISO1(2),ISO1(3),'o', 'MarkerFaceColor','b','MarkerSize',10);
text(ISO1(1)+0.1/sum(ISO1(:)),ISO1(2)+0.1/sum(ISO1(:)),ISO1(3)+0.1/sum(ISO1(:)),'ISO','fontsize',15);
plot3(ISO2(1),ISO2(2),ISO2(3),'o', 'MarkerFaceColor','b','MarkerSize',10);
text(ISO2(1)+0.1/sum(ISO2(:)),ISO2(2)+0.1/sum(ISO2(:)),ISO2(3)+0.1/sum(ISO2(:)),'ISO','fontsize',15);

plot3(C1(1),C1(2),C1(3),'ko', 'MarkerFaceColor','k','MarkerSize',10);
text(C1(1)+0.1/sum(C1(:)),C1(2)+0.1/sum(C1(:)),C1(3)+0.1/sum(C1(:)),'C','fontsize',15);
plot3(C2(1),C2(2),C2(3),'ko', 'MarkerFaceColor','k','MarkerSize',10);
text(C2(1)+0.1/sum(C2(:)),C2(2)+0.1/sum(C2(:)),C2(3)+0.1/sum(C2(:)),'C','fontsize',15);

plot3(LVD1(1),LVD1(2),LVD1(3),'ro', 'MarkerFaceColor','r','MarkerSize',10);
text(LVD1(1)+0.1/sum(LVD1(:)),LVD1(2)+0.1/sum(LVD1(:)),LVD1(3)+0.1/sum(LVD1(:)),'LVD','fontsize',15);
plot3(LVD2(1),LVD2(2),LVD2(3),'ro', 'MarkerFaceColor','r','MarkerSize',10);
text(LVD2(1)+0.1/sum(LVD2(:)),LVD2(2)+0.1/sum(LVD2(:)),LVD2(3)+0.1/sum(LVD2(:)),'LVD','fontsize',15);

plot3(CLVD1(1),CLVD1(2),CLVD1(3),'ro', 'MarkerFaceColor','c','MarkerSize',10);
text(CLVD1(1)+0.1,CLVD1(2)-0.2,CLVD1(3),'CLVD','fontsize',15);
plot3(CLVD2(1),CLVD2(2),CLVD2(3),'ro', 'MarkerFaceColor','c','MarkerSize',10);
text(CLVD2(1)+0.1,CLVD2(2)+0.1,CLVD2(3),'CLVD','fontsize',15);

Two third party functions: sphere3d and arrow3d are from matlab file exchange.

function vec2arc was missing. Now I added here. It is to calculate the arc line defined by two vectors.

function garc=vec2arc(A,B)

theta=acos(dot(A,B));
% normal, binormal, and tangent vectors at point A
N=A; BN=cross(A,B);
T=cross(BN,N); T=T/norm(T);

n=20; t=linspace(0,theta,n)';
N=repmat(N,n,1); T=repmat(T,n,1);
t=repmat(t,1,3);
garc=N.*cos(t)+T.*sin(t);

References:
Hudson, J.A., Pearce, R.G. & Rogers, R.M., 1989. Source time plot for
inversion of the moment tensor, J. geophys. Res., 94(B1), 765–774.
Riedesel, M.A. & Jordan, T.H., 1989. Display and assessment of seismic
moment tensors, Bull. seism. Soc. Am., 79(1), 85–100.
Tape, W. & Tape, C., 2012. A geometric setting for moment tensors,
Geophys. J. Int., in press, doi:10.1111/j.1365-246X.2012.05491.x