This post was motivated by the seismic events starting on 6 January 2020 in Puerto Rico, when the island was devastated by a large (M6.4, 10km depth) earthquake event on the south. A consistent number of strong replicas followed, which left affected citizens wondering whether this was normal.
In this tutorial we’re going to use the Incorporated Research Institutions for Seismology (IRIS) database to download earthquake data to a Matlab variable, where we can work out the statistics of this earthquake sequence.
- Download the IRIS java API here and the IRIS fetch library that allows to connect with the IRIS service.
Change the Matlab working directory to the place where you saved these files. - Add the java API to path:
javaaddpath(‘IRIS-WS-2.0.18.jar’) %change file name accordingly - Define your parameters and ask to retrieve the data on the earthquake events:
Here, we define the minimum magnitude of data we want to get as M=2 and select a rectangular grid in the latitude region of [,] and longitude region of [,], comprising the south of Puerto Rico area. Then, we select our start time as 16 of December 2019, and leave the end time open, so that it can get the most up-to-date data on this database (warning: not real time! do not attempt to use it for earthquake early warning.).
events = irisFetch.Events(‘MinimumMagnitude’,2,…
‘minimumLatitude’,17.5,’maximumLatitude’, 18.5,…
‘minimumLongitude’, -68,’maximumLongitude’, -66,…
‘startTime’,’2019-12-16 00:00:00′) - The results will be saved in the struct events where you can investigate the individual parameters and play around with the data.
Let’s build a plot of these events. - First, just to clean things a bit, we will create independent variables for earthquake times, magnitudes and source depths:
for i=1:length(events)
times(i)=datenum(events(i).PreferredTime);
mags(i)=events(i).PreferredMagnitudeValue;
depths(i)=events(i).PreferredDepth;
end - Now let’s plot the data!
The x-axis will be time and the y-axis will be magnitude, so far so good.But we can also play with the size of each data point so to provide a layman indication of how strongly felt an event is. Such size should increase with magnitude and decrease with depth (shallow earthquakes are felt more strongly). Also some sort of adjustment is needed to the fact that the magnitude scale is logarithmic. So, the simplest choice would be size=exp(mags)/depth.Therefore:
figure;
scatter(times,mags,exp(mags)/depths,’k’) %plot the earthquakes
hold on
scatter(ones(1,length(2:0.01:7))*datenum(datetime(‘now’)),2:0.01:7,’.’,’r’)%draw today’s line
hold off
xlim([min(times) max(times)])
ylim([2 7])
xlabel(‘Time’)
ylabel(‘Magnitude’)
datetick - You should get something like this:
- Now let’s identify the mainshock and check Bath’s law:
(Bath’s law states that the largest aftershock is about 1 magnitude smaller than the mainshock).
%get mainshock parameters:
[~,idx]=max(mags);
disp(strcat(‘Mainshock magnitude M=’, num2str(mags(idx))))
disp(strcat(‘Mainshock time:’, datestr(times(idx))))
%(for local time (AST) subtract datenum(0,0,0,4,0,0))%start to consider something as a true aftershock only after 24h of the mainshock:
aftershock_startTime=times(idx)+datenum(0,0,1,0,0,0);
idxStartTime=find(times<=aftershock_startTime,1);
%get aftershock parameters:
[~,idx_largestaftershock]=max(mags(1:idxStartTime));
disp(strcat(‘Largest aftershock magnitude M=’, …
num2str(mags(idx_largestaftershock))))
disp(strcat(‘Largest aftershock time:’, …
datestr(times(idx_largestaftershock))))
%(for local time (AST) subtract datenum(0,0,0,4,0,0))
%check Bath’s law
if and(mags(idx_largestaftershock)<=1.3,mags(idx_largestaftershock)>=0.8)
disp(‘Baths law is verified, within expected variability’)
else
disp(‘Baths law is not verified for this dataset’)
end - As we’ve seen, Bath’s law is not verified in this particular event.
- We can derive the Gutenberg-Richter relation from the data:
The Gutenberg-Richter law is a relation of the number of earthquakes with the earthquake magnitude, given by N=10a-bM , where N is the number of earthquakes with at least magnitude M, and a and b are constants (10a is the total number of events and b is typically about 1) that would relate to geophysical properties. This is equivalent to state that we should have small earthquakes all the time, but large events should be rare.
i=0;
for M=min(mags):0.1:max(mags)
i=i+1;
N(i)=sum(mags>M);
end
clear M
M=min(mags):0.1:max(mags);
figure
scatter(M,N,’k’)
set(gca,’YScale’,’log’);
xlabel(‘Magnitude’)
ylabel(‘Number of earthquakes with at least this magnitude’)
p=polyfit(M(M<5),log10(N(M<5)),1);
b=p(1);
a=p(2);
disp(strcat(‘Gutenberg-Richter b-value: ‘, num2str(-b)))
disp(strcat(‘Gutenberg-Richter a-value:’, num2str(a))) - You should get something like this:
- The b-value obtained was something around 0.7, which is acceptable for a seismically active area.
As a preliminary conclusion, the Puerto Rico event sequence did follow the typical patterns of a seismically active region. As a suggestion, you may also perform this analysis in a moving spatial grid and obtain a map of the b-value across the island, thus identiying the key seismic regions.
Arthur Valencio 