DAS Earthquake Analysis v2

P. McGill
24Mar25

Investigate an earthquake identified in the Geo-Sense data by Aaron Micallef, as outlined in his 14Mar25 email:

“The 4th February strong seismic event (10.12 am local time on Geo-Sense according to Sintela’s updated timeframe) looks like this on MARS. This is a 2.89 magnitude earthquake at this location 37.1255, -121.525:”


Import a function to read  Sintela DAS V3 format files, which was  developed in the notebook SintelaReadRawV3.nb according to the documentation at https://sintela.atlassian.net/wiki/spaces/DEMO/pages/336917/Phase+Data+Sensor+Data+Sintela+Binary+Format.

sintelaReadRawV3[fileName_String]:= Module[ (* v1.0.0 of this function *)
{fileSizeBytes,inStream,syncWord,hdrSize,hdrVers,sampCnt,sampTime,numChans,numSamps,sampRate,chanSpacing,gaugeLen,trigger,startChan,chanStep,gpsPPSTime,gpsPPSSampOffs,gpsStat,trigFlags,statFlags,packetSizeBytes,totalPackets,totalSamps,packet,header,rdBuf,data={}},

(* Decode header information *)
fileSizeBytes=UnitConvert[FileSize[filename],"Bytes"]//QuantityMagnitude;
inStream=OpenRead[fileName,BinaryFormat->True];
syncWord=BinaryRead[inStream, "UnsignedInteger32"] (* offset 0 *);
hdrSize=BinaryRead[inStream, "UnsignedInteger32"] (* offset 4 *);
hdrVers=BinaryRead[inStream, "UnsignedInteger32"] (* offset 8 *);
If[
syncWord==16^^11223344&&hdrSize==96,
Print["Sintela Data Type 3 - Raw 32 Bit."], 
(Print["Not a Sintela Type 3 File."];Close[inStream];Abort[])
] ;(* end if *)

sampCnt=BinaryRead[inStream, "UnsignedInteger64"] (* offset 12 *);
sampTime=BinaryRead[inStream,"UnsignedInteger64"] (* offset 20 *);
numChans=BinaryRead[inStream,"UnsignedInteger32"] (* offset 28 *);
numSamps=BinaryRead[inStream,"UnsignedInteger32"] (* offset 32 *);
sampRate=BinaryRead[inStream,"Real32"] (* offset 36 *);
chanSpacing=BinaryRead[inStream,"Real32"] (* offset 40 *);
gaugeLen=BinaryRead[inStream,"Real32"] (* offset 44 *);
trigger=BinaryRead[inStream,"UnsignedInteger32"] (* offset 48 *);
startChan=BinaryRead[inStream,"UnsignedInteger32"] (* offset 52 *);
chanStep=BinaryRead[inStream,"UnsignedInteger32"] (* offset 56 *);
gpsPPSTime=BinaryRead[inStream,"UnsignedInteger32"] (* offset 60 *);
gpsPPSSampOffs=BinaryRead[inStream,"Integer32"] (* offset 64 *);
gpsStat=BinaryRead[inStream,"UnsignedInteger32"] (* offset 68 *);
trigFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 72 *);
statFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 73 *);

(* Display first header *)
Print["Packet Channels x Samples : ",numChans," x ",numSamps];Print["Start Channel: ",startChan];
Print["Fs: ",NumberForm[sampRate,{8,3}]," Hz"] ;Print["Channel Spacing: ",NumberForm[chanSpacing,{4,1}]," m"];
Print["Gauge Length: ",NumberForm[gaugeLen,{4,1}]," m"];

Print["Start Time: ",DateString[FromUnixTime[sampTime/10^9,TimeZone->"GMT"]],".",Mod[sampTime,10^9]," GMT"];

packetSizeBytes=(hdrSize+(numChans*numSamps*4));
totalPackets=Round[fileSizeBytes/packetSizeBytes];
totalSamps=totalPackets*numSamps;
Print["Total Samples Per Channel: ",totalSamps];
Print["Total Packets: ",totalPackets];
Print["Bytes per Packet: ",packetSizeBytes];

(* Read and extract data samples (skipping over header blocks) *)
SetStreamPosition[inStream,0];(* Rewind file *)
Monitor[ (* show a progress indicator as we read the packets *)
For[packet=1,packet<=totalPackets,packet++,
(* skip over header and read samples *)
SetStreamPosition[inStream,StreamPosition[inStream]+hdrSize];
rdBuf=BinaryReadList[inStream,"Real32",numChans*numSamps];
AppendTo[data,rdBuf];
] (* end For *),
ProgressIndicator[packet,{1,totalPackets}]
] (* end Monitor *);

(* data are now in sublists of packetSizeBytes length, and we must rearrange to sublists of numChans length *);
data=Partition[Flatten[data],numChans];

Close[inStream];
Return[data]
] (* end module *)

The Geo-Sense DAS had a clock malfunction that caused it to run about 5% fast. To find the earthquake in the DAS data, we have to correct for this drift by using the Excel file “MBARI Geo-sense 1_0 time correction Rev0.03.xlsx” provided by Sintela.

eq1TimeTrue=DateObject[{2025,2,4,18,12,6},TimeZone->"UTC"]

eq1TimeWarped=DateObject[{2025,2,5,3,43,3},TimeZone->"UTC"]

(*filename="/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Geo-Sense\ data/decimator_2025-02-05_03.40.53_UTC_000922.raw";*)

filename="/Volumes/Geo-Sense/2025\ Jan\ Data/decimator_2025-02-05_03.40.53_UTC_000922.raw";

dasDat=sintelaReadRawV3[filename];

file1Strt=DateObject[{2025,2,5,3,40,53},TimeZone->"UTC"]

eq1FileDly=QuantityMagnitude[UnitConvert[DateDifference[file1Strt,eq1TimeWarped],"Seconds"]]

Make a data browser that displays the time near the earthquake, with a 2 Hz highpass filter, for any channel.

Manipulate[ListPlot[HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+25) 1000,chan]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000],Joined->True,PlotRange->{-5,5},DataRange->{0,15},Frame->True,GridLines->Automatic,ImageSize->Large],{chan,1,162,1}]

The earthquake is visible, but the response of the cable varies greatly from channel to channel. A good response is seen on channel 162 at the bitter end of the cable. Plot just that channel with some cursors so we can make some time measurements.

Manipulate[ListPlot[HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+25) 1000,162]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000],Joined->True,PlotRange->{-3,3},DataRange->{0,15},Frame->True,GridLines->Automatic,ImageSize->Large,Prolog->{Thick,Red,Line[{{pWave,-3},{pWave,3}}],Line[{{sWave,-3},{sWave,3}}],Style[Text["Δt = "<>ToString[sWave-pWave]<>" s",{7,1.5}],24]}],{{pWave,3.25},0,7.5,ImageSize->500},{{sWave,10.25},7.5,15,ImageSize->500}]

Using information at

http://eqseis.geosc.psu.edu/cammon/HTML/Classes/IntroQuakes/Notes/waves_and_interior.html

we can find the distance to the earthquake. Note that we had to correct the time delay for the DAS-clock-frequency error.

dt=Quantity[6.3 1.05,"Seconds"]; (* time difference of P- and S-wave arrivals *)
Vp=Quantity[5,"Kilometers"/"Seconds"]; (* P-wave velocity in km/s *)
Vs=Quantity[3,"Kilometers"/"Seconds"]; (* S-wave velocity in km/s *)
dist1=(dt/(((1/Vs)-(1/Vp))))

Find the location of the bitter end of the cable. To do this, import the cable survey data provided by ROV pilot Ethan Brandt.

(*cabSrvy=Drop[#,{3}]&/@Import["/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Cable\ Survey\ Jan\ 2025/GeoSense\ Survey\ Initial.rln","CSV"];*)

cabSrvy=Drop[#,{3}]&/@Import["/Users/mcgill/Documents/My\ Projects/Geo-Sense/Deployment\ Jan\ 2025/Cable\ Survey\ Jan\ 2025/GeoSense\ Survey\ Initial.rln","CSV"];

First@cabSrvy

Length@cabSrvy

Convert the cable survey points from UTM to lat and lon.

cabPts=GeoPosition[GeoGridPosition[#,"UTMZone10"]]&/@cabSrvy;

Take just the last point as the location of the bitter end of the cable.

bitterLoc=Last@cabPts

Find all the earthquakes that occurred within 100 km of that point on 4Feb25.

feb4Eqs=EarthquakeData[bitterLoc,1,{{2025,2,4,0},{2025,2,5,0}},TimeZone->"UTC"]

Get the position of the first earthquake in the list at 1812 UTC.

eq1Loc=First@EarthquakeData[bitterLoc,1,{{2025,2,4,0},{2025,2,5,0}},"Position",TimeZone->"UTC"]["Values"]

Find the distance between the earthquake epicenter and the bitter end of the cable.

UnitConvert[GeoDistance[eq1Loc,bitterLoc],"Kilometers"]

GeoGraphics[{Red,PointSize[Large],Point[eq1Loc],Thick,GeoCircle[eq1Loc,Quantity[#,"Kilometers"]]&/@{2,4,6},Point[bitterLoc],AbsoluteThickness[1],GeoCircle[bitterLoc,dist1]},ImageSize->Full]

Show[Plot3D[((6.3 1.05)/((1/Vs)-(1/Vp))),{Vs,2000,4000},{Vp,4000,6000}],Graphics3D[{Red,PointSize[Large],Point[{3000,5000,49613}]}]]

ArrayPlot[Transpose[BlockMap[RootMeanSquare[#-Mean[#]]&,#,1000]&/@Transpose[dasDat[[100000;;200000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},ColorFunction->"BlueGreenYellow",FrameTicks->Automatic,ImageSize->Full]

ArrayPlot[Transpose[BlockMap[Max[Abs[#-Mean[#]]]&,#,1000]&/@Transpose[dasDat[[100000;;200000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},ColorFunction->"BlueGreenYellow",FrameTicks->Automatic,ImageSize->Full]

To make sure that our excellent results aren’t a fluke, locate a different earthquake that occurred during  the Geo-Sense DAS’s first deployment. List all the events greater than magnitude 2.

EarthquakeData[bitterLoc,2,{{2025,1,28,0},{2025,2,5,23,59,59}},TimeZone->"UTC"]

data=EarthquakeData[bitterLoc,2,{{2025,1,28,0},{2025,2,5,23,59,59}},"Position",TimeZone->"UTC"]["Values"];
GeoGraphics[{Red,PointSize[Large],Point[data],Blue,Point[bitterLoc]}]

data=EarthquakeData[bitterLoc,2,{{2025,1,28,0},{2025,2,5,23,59,59}},"Position",TimeZone->"UTC"]["Values"];
GeoGraphics[{Red,PointSize[Large],Point[data],Blue,Point[bitterLoc]}]

Choose the 1Feb25 00:24:04 event.

eq2TimeTrue=DateObject[{2025,2,1,00,24,4},TimeZone->"UTC"]

eq2TimeWarped=DateObject[{2025,2,1,5,19,59},TimeZone->"UTC"]

(*filename2="/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Geo-Sense\ data/decimator_2025-02-01_05.15.23_UTC_000481.raw";*)

filename2="/Volumes/Geo-Sense/2025\ Jan\ Data/decimator_2025-02-01_05.15.23_UTC_000481.raw";

dasDat2=sintelaReadRawV3[filename2];

file2Strt=DateObject[{2025,2,1,5,15,23},TimeZone->"UTC"]

eq2FileDly=QuantityMagnitude[UnitConvert[DateDifference[file2Strt,eq2TimeWarped],"Seconds"]]

Make a data browser that displays the time near the earthquake, with a 2 Hz highpass filter, for any channel.

Manipulate[ListPlot[HighpassFilter[dasDat2[[(eq2FileDly+5) 1000;;(eq2FileDly+20) 1000,chan]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000],Joined->True,PlotRange->{-1,1},DataRange->{0,15},Frame->True,GridLines->Automatic,ImageSize->Large],{chan,1,162,1}]

Manipulate[ListPlot[dasDat2[[(eq2FileDly+50) 1000;;(eq2FileDly+75) 1000,chan]],Joined->True,PlotRange->{-5,5},DataRange->{0,25},Frame->True,GridLines->Automatic,ImageSize->Large],{chan,1,162,1}]

Manipulate[ListPlot[HighpassFilter[dasDat2[[(eq2FileDly+5) 1000;;(eq2FileDly+15) 1000,162]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000],Joined->True,PlotRange->{-1,1},DataRange->{0,10},Frame->True,GridLines->Automatic,ImageSize->Large,Prolog->{Thick,Red,Line[{{pWave,-3},{pWave,3}}],Line[{{sWave,-3},{sWave,3}}],Style[Text["Δt = "<>ToString[sWave-pWave]<>" s",{5,0.5}],24]}],{{pWave,2.5},0,5,ImageSize->500},{{sWave,7.5},5,10,ImageSize->500}]

dt=Quantity[4.12  1.051,"Seconds"]; (* time difference of P- and S-wave arrivals *)
Vp=Quantity[5,"Kilometers"/"Seconds"]; (* P-wave velocity in km/s *)
Vs=Quantity[3,"Kilometers"/"Seconds"]; (* S-wave velocity in km/s *)
dist2=(dt/(((1/Vs)-(1/Vp))))

eq2Loc=EarthquakeData[(Entity["ncnc75126526", Entity["Earthquake","ncnc75126526"], "Entity[\"Earthquake\", \"ncnc75126526\"]", "Earthquake"]),"Position",TimeZone->"UTC"]

Find the distance between the earthquake epicenter and the bitter end of the cable.

UnitConvert[GeoDistance[eq2Loc,bitterLoc],"Kilometers"]

GeoGraphics[{Red,PointSize[Large],Point[eq2Loc],Thick,GeoCircle[eq2Loc,Quantity[#,"Kilometers"]]&/@{1,2,3},Point[bitterLoc],AbsoluteThickness[1],GeoCircle[bitterLoc,dist2]},ImageSize->Full]

Import the GCC-PHAT TDOA function developed in the notebook “DAS GCC-PHAT.nb”.

gccPhat[sig_List,ref_List]:=Module[{sigt,reft,prod,denom,delay}, (* version 1.0.0 of function *)
If[Length@sig!=Length@ref,Print["Error: Signal and reference must be same length"];Return[]];
sigt=Fourier[Join[sig,Table[0,Length[sig]]]];(* transform to frequency domain *)
reft=Fourier[Join[ref,Table[0,Length[ref]]]];
prod=sigt Conjugate[reft]; (* multiply transforms *)
denom=Clip[Abs[prod],{(10^(-6)),Infinity}]; (* find normalization and clip small values *)
delay=Abs[InverseFourier[prod/denom]];
Return[Join[delay[[Length[delay]/2+2;;]],delay[[;;Length[delay]/2]]]]
] (* end Module *)

According to Aaron and Rebecca’s cable survey map, the cable is pretty straight from about 100 to 300 m away from the DAS (eastern) end of the cable.

Round[(({100,300})/6.3)]

This is about channel 16 to channel 48.

eq1chan16=HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+50) 1000,16]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000];
ListLinePlot[eq1chan16,PlotRange->{-5,5},DataRange->{0,40},Frame->True,GridLines->Automatic,ImageSize->Large]

eq1chan48=HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+50) 1000,48]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000];
ListLinePlot[eq1chan48,PlotRange->{-5,5},DataRange->{0,40},Frame->True,GridLines->Automatic,ImageSize->Large]

Length@eq1chan16

eq1xcor1=gccPhat[eq1chan48[[;;-2]],eq1chan16[[;;-2]]];
ListLinePlot[eq1xcor1[[39950;;40050]],PlotRange->All,DataRange->{-50,50}]

This shows no delay. Find what the delay would be between two channels 32 channels apart for energy travelling at 5 km/s in an end-fire orientation.

((32 6.3)/5000)

We should see a 30- to 40-sample delay but it’s not evident in the original data or the cross correlation.

Let’s find the delay between the earliest channel where the earthquake is visible, channel 4, and the bitter end of the cable, channel 162.

eq1chan4=HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+50) 1000,4]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000];
ListLinePlot[eq1chan4,PlotRange->{-5,5},DataRange->{0,40},Frame->True,GridLines->Automatic,ImageSize->Large]

eq1chan162=HighpassFilter[dasDat[[(eq1FileDly+10) 1000;;(eq1FileDly+50) 1000,162]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000];
ListLinePlot[eq1chan162,PlotRange->{-5,5},DataRange->{0,40},Frame->True,GridLines->Automatic,ImageSize->Large]

eq1xcor2=gccPhat[eq1chan162[[;;-2]],eq1chan4[[;;-2]]];
ListLinePlot[eq1xcor2[[39950;;40050]],PlotRange->All,DataRange->{-50,50}]

There’s no apparent delay between channels 4 and 162. Is this a flaw in the GCC-PHAT code? Try finding the delay between channel 162 and an artificially delayed version of channel 162.

dlyeq1chan162=Drop[Join[Table[0,25],eq1chan162],-25];
ListLinePlot[{eq1chan162,dlyeq1chan162},PlotRange->{{3800,4800},All}]

eq1xcor3=gccPhat[dlyeq1chan162[[;;-2]],eq1chan162[[;;-2]]];
ListLinePlot[eq1xcor3[[39950;;40050]],PlotRange->All,DataRange->{-50,50}]

The following are experiments in using the waveform envelope to detect the onset of the s-wave.

ListPlot[HighpassFilter[dasDat2[[(eq2FileDly+5) 1000;;(eq2FileDly+15) 1000,162]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000],Joined->True,PlotRange->{-1,1},DataRange->{0,10},Frame->True,GridLines->Automatic,ImageSize->Large]

seisDat=HighpassFilter[dasDat2[[(eq2FileDly+5) 1000;;(eq2FileDly+15) 1000,162]],(Quantity[2, "Hz", hertz, "Hertz"]),1001,SampleRate->1000];

ListPlot[{seisDat,Abs[seisDat+I 3 HilbertFilter[seisDat,(Quantity[7, "Hz", hertz, "Hertz"]),SampleRate->1000]]},DataRange->{0,10},PlotRange->{-1,1},Joined->True]

ListLinePlot[BandpassFilter[dasDat2[[(eq2FileDly+5) 1000;;(eq2FileDly+15) 1000,162]],{(Quantity[2, "Hz", hertz, "Hertz"]),(Quantity[10, "Hz", hertz, "Hertz"])},1001,SampleRate->1000],Frame->True,GridLines->Automatic,ImageSize->Large]

sta=MovingAverage[seisDat,100];

lta=MovingAverage[seisDat,1000];

ListLinePlot[sta[[450;;9902-451]]/lta,Frame->True,GridLines->Automatic,ImageSize->Large]

Length@sta

Length@lta

Length@seisDat

Length@sta[[450;;9902-451]]