This note documents a simple in-class activity that I use in my GEOSC 109H earthquakes and society class. In class I describe the Gutenberg-Richter Relation that represents the distribution of earthquakes as a function of magnitude (a nice paper-based activity is available from SCEC). I show the students a plot from the GCMT catalog in class, in this activity, they use Mathematica to download data for a few randomly chosen years and see if the pattern I showed in class holds up.

The Mathematica notebook begins with a few utility functions that allow students to use web services to acquire seismicity data and a few functions to list some events (which is partly a check to see that things are working). Below is a listing of the search script (an example usage is provided below).

comcatRectangularSearchPath[start_, end_, minlat_, maxlat_, minlon_, maxlon_] := Module[{},
  "http://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=" <> 
   DateString[start, "ISODateTime"] <> "&endtime=" <> 
   DateString[end, "ISODateTime"] <> 
   "&minmagnitude=4.0&minlatitude=" <> 
   ToString[Round[minlat, 0.001]] <>
   "&maxlatitude=" <> 
   ToString[Round[maxlat, 0.001]] <>
   "&minlongitude=" <> 
   ToString[Round[minlon, 0.001]] <> "&maxlongitude=" <> 
   ToString[Round[maxlon, 0.001]] <> "&orderby=time-asc"
  ]

comcatRectangularSearch[start_, end_, minlat_: - 90, maxlat_: 90, 
  minlon_: - 180, maxlon_: 180] := Module[{path},
   path = comcatRectangularSearchPath[start, end, minlat, maxlat, minlon, maxlon];
  {path, Import[path, "CSV"]}
  ]

The scripts are simplistic, and to avoid errors for some of the more active years, I limit the magnitude range of the analysis to events greater than or equal to magnitude 4.0. The code below is a function to list the largest events in the search results.

listEventHeadings = {"Origin Time", "Latitude", "Longitude", "Depth", 
   "Magnitude", "MagType"};
(**)
EventSummaryString[eq_] := Module[{s0, s1},
  s0 = DateString[eq[[1]]] <> "\t" 
    <> ToString[eq[[2]]] <> "\t"
    <> ToString[eq[[3]]] <> "\t"
    <> ToString[eq[[4]]] <> "\t"
    <> ToString[eq[[5]]] <> "\t"
    <> ToString[eq[[6]]] <> "\t"
  ]
(**)
ListNLargestEvents[cat_, n_] := Module[{sorted, nevents},
  sorted = Sort[cat, #1[[5]] > #2[[5]] &];
  nevents = {#[[1]], #[[2]], #[[2]], #[[4]], #[[5]], #[[6]]} & /@ 
    sorted[[1 ;; n]];
  TableForm[nevents, TableHeadings -> {None, listEventHeadings}]
  ]

Using the Scripts

Students perform searches on individual years using the following commands. The search function is called and it returns a two-item list containing the search string (an http address useful for debugging and the output of the web services search). For convenience, I strip out the origin time and magnitude values of the events and store them in the lists called otimes and mags.

theYEAR = 2005;
(**)
usgsresponse = comcatRectangularSearch[{theYEAR, 1, 1}, {theYEAR + 1, 1, 1, 0, 0}];
headings = usgsresponse[[2]][[1]]
usgscat = Drop[usgsresponse[[2]], 1]

otimes = #[[1]] & /@ usgscat;
mags = #[[5]] & /@ usgscat;

To perform a search for a different year, the just have to edit the cell to change the value of theYEAR. To see a list of the five largest events, the next cell is

(* List the five largest events of the year *)
ListNLargestEvents[usgscat, 5]

The output of that call in this case is (the formatting is better in the notebook, things line up).

Origin Time	Latitude	Longitude	Depth	Magnitude	MagType
2005-03-28T16:09:36.530Z	2.085	2.085	30	8.6	mww
2005-06-13T22:44:33.900Z	-19.987	-19.987	115.6	7.8	mwb
2005-10-08T03:50:40.800Z	34.539	34.539	26	7.6	mwc
2005-09-09T07:26:43.730Z	-4.539	-4.539	90	7.6	mwc
2005-09-26T01:55:37.670Z	-5.678	-5.678	115	7.5	mwb

The Gutenberg-Richter Analysis

I perform the G-R analysis using the classic table listing the numbers of events with magnitude greater than or equal to a sequence of magnitudes. The cumulative distribution is probably a better way to construct the model, but this is a 100-level class for non-science students, the intuitive version is preferred. The code below will tally the numbers that we need, starting a M 4.0 and using bin widths of 0.5 magnitude units.

(* Compute the Gutenberg Richter Statistics *)
m = 4.0;
dm = 0.5;
grStats = {};
While[m <= 8.5,
 ngtm = Length[Select[mags, # >= m &]];
 grStats = Append[grStats, {m, ngtm}];
 m += dm;
 ]
grStats;
counts = Transpose[grStats][[2]];
bcenters = Transpose[grStats][[1]];
(**)
Print["Global earthquake Statistics for the year: " <> ToString[theYEAR] <>"."]
TableForm[{Range[1, Length[counts]], bcenters, counts} // Transpose, 
  TableHeadings -> {None, {"Index", "Magnitude","N \[GreaterEqual] Magnitude"}}, 
  TableAlignments -> Center] // Print

Global earthquake Statistics for the year: 2005.
Index	Magnitude	N >= Magnitude
1	    4.	           15771
2	    4.5	           7886
3	    5.	           1843
4	    5.5	           533
5	    6.	           151
6	    6.5	           45
7	    7.	           11
8	    7.5	           5
9	    8.	           1
10	    8.5	           1

I list the index because to choose the range over which to fit the line to the pattern, we need to insure that none of the bins we use are empty. By default, I start with the second bin and include bins 2 through 8. This seems to work well, but to be honest, having something not work is a good lesson in computing for the students. Next I create a plot of the values, N vs Magnitude.

(**)
mark = Graphics[{EdgeForm[GrayLevel[0.05]], Opacity[0.65], FaceForm[Red], 
    Disk[{0, 0}, 0.05], PlotRangeClipping -> False}];
(**)
pplt = ListPlot[{bcenters, Log10[counts]} // Transpose,
   Frame -> True, Axes -> False,
   FrameLabel -> {"Magnitude", 
     "\!\(\*SubscriptBox[\(Log\), \(10\)]\)[N \[GreaterEqual] M]"},
   RotateLabel -> True, ImageSize -> 512,
   BaseStyle -> {24, FontFamily -> "Helvetica"}, 
   PlotMarkers -> {mark, 0.04}, PlotRange -> {{3, 9.}, {-1, 5}}];
(**)
pplt // Print

Here’s the output

Next we fit a line to the values. I used built-in commands a simply write the slope and intercept values to the notebook. Students can copy and paste the information into an online activity form (in the course management system we use, Canvas).

i0 = 2;
i1 = 8;
Clear[x];
lm = LinearModelFit[{bcenters[[i0 ;; i1]], 
     Log10[counts[[i0 ;; i1]]]} // Transpose, x, x];
(**)
p = lm["BestFitParameters"];
(**)
Print[Style["The best fit line slope: " <> ToString[p[[2]]] <> 
   " intercept: " <> ToString[p[[1]]], {18, Bold}]]

Here’s the output

The best fit line slope: -1.07967 intercept: 8.68685

We can overlay the line onto the data plot using

(**)
lplt = Plot[lm[x], {x, 3, 9}, 
       PlotStyle -> {RGBColor[0.5, 0.5, 0.5], Thickness[0.01]}];
(**)
grplt = Show[pplt, lplt,
   FrameStyle -> AbsoluteThickness[0.85],
   ImageSize -> 512, 
   PlotLabel ->  ToString["Global EQ Activity, Year: " <> ToString[theYEAR]],
   GridLines -> Automatic,
   Epilog -> {Text[ Style["Data Source\nUSGS", 18, Background -> White], 
    {8,3.5}, {1, -1}]}
   ];
(**)
grplt // Print

The display is more formal and summarizes the analysis.

A Global Earthquake Map

I finish with a simple map of the earthquakes. I don’t show the background because I want students to focus on the spatial patterns in earthquake locations. Repeating this for a few years communicates the ideas that the first-order earthquake patterns are the same each year, but each year is different.

elocs = {#[[2]], #[[3]], #[[5]]} & /@ Select[usgscat, #[[5]] > 4 &];

west = -180;
east = 180;
south = -90;
north = 90;

eqradius[m_] := N[10^(0.5*m - 2.25)];
scaleFactor = 10;

GeoGraphics[{GeoStyling[Opacity[0.5], FaceForm[{Red}],EdgeForm[{Black, Thin}]], 
   GeoDisk[{#[[1]], #[[2]]}, 
   Quantity[scaleFactor*Max[{eqradius[#[[3]]], 4}]  , "km"]] & /@ elocs},
  GeoRange -> {{south, north}, {west, east}},
  BaseStyle -> {21, FontFamily -> "Helvetica"},
  GeoProjection -> {"Robinson", "Centering" -> 0, "GridOrigin" -> {0, 0}},
  FrameTicks -> {{None, None}, {Automatic, Automatic}},
  GeoBackground -> None,
  GeoGridLines -> {Range[south, north, 15.0], Range[west, east, 15.0]},
  GeoGridLinesStyle -> {{Opacity[.25], Black}, {Opacity[.25], Black}},
  ImageSize -> 800] // Print

Here’s is the result for 2005. Symbol size is proportional to rupture size, but I exaggerated the the size by a factor of three to make smaller magnitude events visible. I also limit the minimum earthquake radius to 4 km, so the small events are visible.

Note that Mathematica plots circular regions onto the map projection – this leads to an important discussion of map projections and their distortions.

Summary

Mathematica provides useful tools for searching, analyzing and displaying earthquake locations. I used the above code to have students check the validity of the Gutenberg-Richter Relation for earthquakes, stressing the evidence-based nature of this important empirical pattern. We return to similar analyses for individual large earthquakes later in the semester.