Geologist & seismologist: MQL recipes

Recipes

P and S arrivals and epicentral distance
Energy of a magnitude 9.1 earthquake
b-value of an earthquake catalogue
Radiocarbon age with its uncertainty
Rb-Sr isochron age of a granite
Arrival delay between two seismic stations
Depth, mass and size of an ore body from a gravity anomaly
Palaeolatitude and virtual geomagnetic pole of a lava flow
Thermal thickness and heat flow of the oceanic lithosphere
Decomposition of an earthquake moment tensor

executed on MentDB

P and S arrivals and epicentral distance

Band-pass a synthetic seismogram to remove the microseism, pick the P and S arrivals with a peak detector and get the distance from the S-P time.

Pipeline

01 · recordSeismogram at 50 Hznum table→02 · filterButterworth band-passdsp filter→03 · pickMain peaksdsp peaks→04 · distanceS-P time to kmmath calc

Script

admin

#1. 20 s seismogram at 50 Hz: microseism 0.2 Hz, P wave at 6 Hz, S wave at 3 Hz;
-> "[sig]" "0.8*sin(2*pi*0.2*t)+1.0*exp(-((t-5.6)/0.4)^2)*sin(2*pi*6*t)+2.5*exp(-((t-12.3)/0.7)^2)*sin(2*pi*3*t)";
-> "[x]" (matrix col (num table [sig] "t" 0 20 0.02) 2);
#2. Band-pass 1 to 10 Hz Butterworth to remove the microseism;
-> "[f]" (dsp filter [x] 50 "bandpass" "[1, 10]" 4);
#3. Main peaks at least 3 s apart;
-> "[pk]" (dsp peaks [f] "{\"height\": 0.5, \"distance\": 150}");
json load "rx_seis" [pk];
-> "[tP]" (math calc (concat (json select "rx_seis" /indices[0]) "/50"));
-> "[tS]" (math calc (concat (json select "rx_seis" /indices[1]) "/50"));
#4. S-P time and epicentral distance with vp 6.0 and vs 3.5 km/s;
concat "P at " [tP] " s, S at " [tS] " s, S-P = " (math round (math calc "[tS]-[tP]") 2) " s, distance = " (math round (math calc "([tS]-[tP])*6.0*3.5/(6.0-3.5)") 1) " km"

Result

mentdb

P at 5.54 s, S at 12.42 s, S-P = 6.88 s, distance = 57.8 km

Go further

Use three stations and num solve_system to locate the epicentre
Look at the spectrum of each phase with dsp spectrum

executed on MentDB

Energy of a magnitude 9.1 earthquake

Compute the radiated energy and the seismic moment of the Tohoku earthquake, its TNT equivalent and how many M6 shocks it is worth.

Pipeline

01 · energylog10 E = 1.5 M + 4.8math calc→02 · convertJ to MWhunit convert→03 · compareMoment and ratiomath calc

Script

admin

#1. Moment magnitudes of the Tohoku 2011 earthquake and of a moderate shock;
-> "[M1]" 9.1;
-> "[M2]" 6.0;
#2. Radiated energy, Gutenberg-Richter log10 E = 1.5 M + 4.8 in joules;
-> "[E1]" (math calc "10^(1.5*[M1]+4.8)");
-> "[E2]" (math calc "10^(1.5*[M2]+4.8)");
#3. Seismic moment in N.m, TNT equivalent and energy ratio;
-> "[M0]" (math calc "10^(1.5*[M1]+9.1)");
concat "E = " (math round (math calc "[E1]/1e18") 2) "e18 J = " (math round (math calc "[E1]/4.184e15") 0) " Mt TNT = " (math round (math calc (concat (unit convert [E1] "J" "MWh") "/1e6")) 0) " TWh, M0 = " (math round (math calc "[M0]/1e22") 1) "e22 N.m, ratio to M" [M2] " = " (math round (math calc "[E1]/[E2]") 0)

Result

mentdb

E = 2.82e18 J = 674 Mt TNT = 783 TWh, M0 = 5.6e22 N.m, ratio to M6.0 = 44668

Go further

Get the moment magnitude from a measured seismic moment
Tabulate the energy for M 4 to 9 with num table

executed on MentDB

b-value of an earthquake catalogue

Fit the Gutenberg-Richter law to the cumulative counts of a regional catalogue, then estimate the yearly rate and the return period of M >= 6 events.

Pipeline

01 · fitExponential fit of Nstat fit→02 · b-valueb and a from the fitmath log10→03 · hazardRate and return periodmath calc

Script

admin

#1. Regional catalogue: number of earthquakes with magnitude >= M;
-> "[M]" "[2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]";
-> "[N]" "[1420, 468, 151, 52, 17, 5, 2]";
#2. Gutenberg-Richter law log10 N = a - b M as an exponential fit;
-> "[fit]" (stat fit [M] [N] "exp");
json load "rx_gr" [fit];
-> "[b]" (math calc (concat "-1*" (json select "rx_gr" /b) "/ln(10)"));
-> "[a]" (math log10 (json select "rx_gr" /a));
#3. Expected yearly rate of M >= 6 events and return period, catalogue of 10 years;
-> "[n6]" (math calc "10^([a]-[b]*6)/10");
concat "b-value = " (math round [b] 3) ", a = " (math round [a] 3) ", M>=6 rate = " (math round [n6] 4) " per year, return period = " (math round (math calc "1/[n6]") 0) " years"

Result

mentdb

b-value = 0.96, a = 5.068, M>=6 rate = 0.0202 per year, return period = 49 years

Go further

Use the maximum likelihood estimate of Aki on the raw magnitudes
Get the Poisson probability of at least one event in 50 years

executed on MentDB

Radiocarbon age with its uncertainty

From the fraction of modern carbon of a charcoal sample, compute the conventional radiocarbon age and the age with the true half-life, with propagated uncertainties.

Pipeline

01 · ageLibby age -8033 ln Fnum uncertainty→02 · half-lifeAge with T1/2 = 5730 +- 40num uncertainty→03 · readText and error budgetjson select

Script

admin

#1. Charcoal sample: fraction of modern carbon F = 0.6123 +- 0.0035;
-> "[F]" "[0.6123, 0.0035]";
#2. Conventional radiocarbon age with the Libby mean life 8033 years, with its uncertainty;
-> "[age]" (num uncertainty "-8033*ln(F)" (concat "{\"F\": " [F] "}"));
json load "rx_c14" [age];
#3. Same fraction with the true half-life 5730 +- 40 years;
-> "[true]" (num uncertainty "-T/ln(2)*ln(F)" (concat "{\"F\": " [F] ", \"T\": [5730, 40]}"));
json load "rx_c14t" [true];
concat "radiocarbon age = " (json select "rx_c14" /text) " BP, with the 5730 yr half-life = " (json select "rx_c14t" /text) " yr, share of the half-life in the error = " (math round (math calc (concat "100*" (json select "rx_c14t" /variance_share/T))) 1) " %"

Result

mentdb

radiocarbon age = 3940 ± 50 BP, with the 5730 yr half-life = 4060 ± 60 yr, share of the half-life in the error = 26.4 %

Go further

Calibrate the age on a curve with num interpolate
Date a potassium-argon sample the same way

executed on MentDB

Rb-Sr isochron age of a granite

Fit the isochron line through the whole-rock ratios, get the initial strontium ratio and the slope with its error, and turn it into an age in Ma.

Pipeline

01 · fitLine with parameter errorsstat fit_model→02 · readSlope and errorjson select→03 · ageln(1+s)/lambdanum uncertainty

Script

admin

#1. Whole-rock samples of a granite: 87Rb/86Sr and 87Sr/86Sr ratios;
-> "[x]" "[0.25, 0.80, 1.45, 2.10, 3.05]";
-> "[y]" "[0.70862, 0.71841, 0.72931, 0.74082, 0.75668]";
#2. Isochron line y = a + s*x with the errors of the parameters;
-> "[fit]" (stat fit_model "a+s*x" "x" [x] [y] "{\"a\": 0.7, \"s\": 0.01}");
json load "rx_iso" [fit];
-> "[s]" (json select "rx_iso" /parameters/s);
-> "[us]" (json select "rx_iso" /errors/s);
#3. Age t = ln 1+s / lambda with lambda 1.42e-11 per year, and its uncertainty;
-> "[t]" (num uncertainty "ln(1+s)/1.42e-11/1e6" (concat "{\"s\": [" [s] ", " [us] "]}"));
json load "rx_isot" [t];
concat "initial 87Sr/86Sr = " (math round (json select "rx_iso" /parameters/a) 5) ", slope = " (math round [s] 6) ", age = " (json select "rx_isot" /text) " Ma"

Result

mentdb

initial 87Sr/86Sr = 0.70451, slope = 0.017163, age = 1198 ± 8 Ma

Go further

Weight the points by their analytical errors
Do the same for a Sm-Nd isochron

executed on MentDB

Arrival delay between two seismic stations

The same P wavelet is recorded by two stations 12 km apart. The cross-correlation of the two records peaks at the delay between them, which gives the apparent velocity of the wave across the pair.

Pipeline

01 · recordsTwo records at 50 Hznum table→02 · xcorrCross-correlationdsp correlate→03 · lagLag of the maximumdsp peaks→04 · velocityApparent velocitymath calc

Script

admin

#1. Same P wavelet recorded by two stations 12 km apart, 50 samples/s over 10 s;
-> "[A]" (matrix col (num table "exp(-((t-3)/0.15)^2)*sin(2*pi*5*t)" "t" 0 10 0.02) 2);
-> "[B]" (matrix col (num table "0.6*exp(-((t-5)/0.15)^2)*sin(2*pi*5*(t-2))" "t" 0 10 0.02) 2);
#2. Cross-correlation B against A, lags from -500 to 500 samples;
-> "[c]" (dsp correlate [B] [A]);
#3. Lag of the maximum: index - len of A - 1;
json load "rx_xc" (dsp peaks [c] "{\"distance\": 2000}");
-> "[lag]" (math calc (concat "(" (json select "rx_xc" /indices[0]) "-500)*0.02"));
#4. Apparent velocity across the pair of stations;
concat "delay = " (math round [lag] 3) " s, apparent velocity = " (math round (math calc "12/[lag]") 2) " km/s"

Result

mentdb

delay = 2 s, apparent velocity = 6 km/s

Go further

With three stations the two delays give the azimuth of the source (beam forming).
Filter the records first with dsp filter bandpass when they are noisy.

executed on MentDB

Depth, mass and size of an ore body from a gravity anomaly

A gravity profile over a buried ore body shows a peak and a half-width. With the sphere model, the script gets the depth of the centre, the excess mass, the radius with a cube root for a known density contrast, and checks the anomaly at the half-width with the forward formula.

Pipeline

01 · constantG of CODATAconst get→02 · inverseDepth and excess massmath calc→03 · radiusCube rootmath cbrt→04 · forwardDistance to the centremath hypot

Script

admin

#1. Gravity survey over a buried ore body: peak anomaly 0.42 mGal, half-width of the profile 460 m;
-> "[gmax]" (math calc "0.42e-5"); -> "[xh]" 460;
-> "[G]" (const get "G"); -> "[drho]" 1500;
#2. Sphere model: depth of the centre z = x1/2 / 0.766 and excess mass M = gmax z^2 / G;
-> "[z]" (math calc "[xh]/0.766");
-> "[M]" (math calc "[gmax]*[z]^2/[G]");
#3. Radius of a sphere with a density contrast of 1500 kg/m3, cube root;
-> "[R]" (math cbrt (math calc "3*[M]/(4*pi*[drho])"));
#4. Forward check at x = 460 m: g = G M z / distance^3;
-> "[d]" (math hypot [xh] [z]);
-> "[g]" (math calc "[G]*[M]*[z]/[d]^3");
concat "depth = " (math round [z] 0) " m, excess mass = " (math round_sig (math calc "[M]/1e9") 3) " Mt, radius = " (math round [R] 1) " m, anomaly at the half-width = " (math round (math calc "[g]*1e5") 3) " mGal, ratio = " (math round (math calc "[g]/[gmax]") 3) ", top of the body at " (math round (math calc "[z]-[R]") 0) " m"

Result

mentdb

depth = 601 m, excess mass = 22.7 Mt, radius = 153.4 m, anomaly at the half-width = 0.21 mGal, ratio = 0.5, top of the body at 447 m

Go further

Fit the whole profile with stat fit_model on G M z / (x^2 + z^2)^1.5 to get z and M with their errors.
Compare with a horizontal cylinder: half-width = z, anomaly in 1/(x^2+z^2).

executed on MentDB

Palaeolatitude and virtual geomagnetic pole of a lava flow

From the mean palaeomagnetic declination and inclination of a lava flow, the script applies the dipole formula to get the palaeolatitude, the magnetic colatitude with the cotangent, the position of the virtual geomagnetic pole and the intensity of the dipole field at that latitude.

Pipeline

01 · dipoletan I = 2 tan lambdamath atan→02 · colatitudep with the cotangentmath cot→03 · poleVGP latitude and longitudemath calc→04 · anglesDegreesmath rad_to_deg

Script

admin

#1. Permian lava flow sampled at 45.2 N, 3.1 E: mean palaeomagnetic direction D = 348 deg, I = 38 deg;
-> "[ls]" (math deg_to_rad 45.2); -> "[phis]" (math deg_to_rad 3.1);
-> "[D]" (math deg_to_rad 348); -> "[I]" (math deg_to_rad 38);
#2. Dipole formula tan I = 2 tan lambda: palaeolatitude, and magnetic colatitude p with the cotangent;
-> "[lam]" (math atan (math calc (concat (math tan [I]) "/2")));
-> "[p]" (math atan (math calc (concat "2*" (math cot [I]))));
#3. Virtual geomagnetic pole: latitude, then longitude difference;
-> "[lp]" (math calc "asin(sin([ls])*cos([p])+cos([ls])*sin([p])*cos([D]))");
-> "[beta]" (math calc "asin(sin([p])*sin([D])/cos([lp]))");
#4. Dipole field intensity at the palaeolatitude, B0 = 30 microtesla at the equator;
-> "[B]" (math calc "30*sqrt(1+3*sin([lam])^2)");
concat "palaeolatitude = " (math round (math rad_to_deg [lam]) 2) " deg, colatitude p = " (math round (math rad_to_deg [p]) 2) " deg, VGP = " (math round (math rad_to_deg [lp]) 2) " N, " (math round (math calc (concat "mod(" (math rad_to_deg (math calc "[phis]+[beta]")) "+360, 360)")) 2) " E, dipole field = " (math round [B] 2) " microtesla"

Result

mentdb

palaeolatitude = 21.34 deg, colatitude p = 68.66 deg, VGP = 64.18 N, 336.7 E, dipole field = 35.46 microtesla

Go further

Add the alpha95 cone of confidence to get dp and dm of the pole.
Compare the pole with the apparent polar wander path of Europe to date the flow.

executed on MentDB

Thermal thickness and heat flow of the oceanic lithosphere

With the half-space cooling model of a 50 Myr old oceanic plate, the script gets the temperature at 30 km depth with the error function, the thickness of the thermal lithosphere with the inverse error function, the surface heat flow and the ocean depth of the Parsons-Sclater law.

Pipeline

01 · ageMyr to secondsunit convert→02 · profileT = Tm erfmath erf→03 · thicknessDepth of 0.9 Tmmath erf_inv→04 · flowHeat flow and depthmath calc

Script

admin

#1. Oceanic lithosphere 50 Myr old: mantle at 1350 C, diffusivity 1e-6 m2/s, conductivity 3.3 W/m/K;
-> "[Tm]" 1350; -> "[kappa]" 1e-6; -> "[k]" 3.3;
-> "[t]" (unit convert 50e6 "yr" "s");
-> "[L]" (math calc "2*sqrt([kappa]*[t])");
#2. Half-space cooling: temperature at 30 km depth, T = Tm erf z / 2 sqrt kappa t;
-> "[T30]" (math calc (concat "[Tm]*" (math erf (math calc "30e3/[L]"))));
#3. Thickness of the thermal lithosphere, the depth where T = 0.9 Tm, with the inverse error function;
-> "[zL]" (math calc (concat "[L]*" (math erf_inv 0.9)));
#4. Surface heat flow and ocean depth of the Parsons-Sclater law 2500 + 350 sqrt t;
-> "[q]" (math calc "[k]*[Tm]/sqrt(pi*[kappa]*[t])");
concat "T at 30 km = " (math round [T30] 0) " C, lithosphere thickness = " (math round (math calc "[zL]/1000") 1) " km, heat flow = " (math round (math calc "[q]*1000") 1) " mW/m2, ocean depth = " (math round (math calc "2500+350*sqrt(50)") 0) " m"

Result

mentdb

T at 30 km = 549 C, lithosphere thickness = 92.4 km, heat flow = 63.3 mW/m2, ocean depth = 4975 m

Go further

Loop on the age from 0 to 150 Myr and plot the depth to see where the plate model flattens.
Compare the heat flow with the measured values near a ridge, lowered by hydrothermal circulation.

executed on MentDB

Decomposition of an earthquake moment tensor

A centroid moment tensor is split into an isotropic part and a deviatoric part. The script gets the principal T, N and P axes as eigenvalues, the double-couple percentage, the scalar moment from the Frobenius norm and the moment magnitude.

Pipeline

01 · isotropicTrace / 3matrix trace→02 · deviatoricM - iso Imatrix diag→03 · axesEigenvaluesmatrix eigen→04 · momentScalar momentmatrix norm_frobenius→05 · magnitudeMwmath calc

Script

admin

#1. Centroid moment tensor of an earthquake, units 1e17 N m, r-theta-phi frame;
-> "[M]" "[[-1.2, 0.4, 2.9], [0.4, 3.5, -1.1], [2.9, -1.1, -1.6]]";
#2. Isotropic part: trace / 3 on a diagonal matrix, deviatoric part;
-> "[iso]" (math calc (concat (matrix trace [M]) "/3"));
-> "[dev]" (matrix sub [M] (matrix diag (concat "[" [iso] "," [iso] "," [iso] "]")));
#3. Principal axes: eigenvalues T, N, P of the deviatoric tensor;
json load "rx_p4e_eig" (matrix eigen [dev]);
-> "[l1]" (json select "rx_p4e_eig" "/values/[0]");
-> "[l2]" (json select "rx_p4e_eig" "/values/[1]");
-> "[l3]" (json select "rx_p4e_eig" "/values/[2]");
json unload "rx_p4e_eig";
#4. Double-couple percentage 1 - 2 eps, scalar moment from the Frobenius norm, moment magnitude;
-> "[eps]" (math calc "abs([l2])/max(abs([l1]), abs([l3]))");
-> "[M0]" (math calc (concat (matrix norm_frobenius [dev]) "/sqrt(2)*1e17"));
-> "[Mw]" (math calc "2/3*(log10([M0])-9.1)");
concat "isotropic = " (math round [iso] 4) ", T N P axes = " (math round [l1] 3) ", " (math round [l2] 3) ", " (math round [l3] 3) ", double couple = " (math round (math calc "100*(1-2*[eps])") 1) " %, M0 = " (math round_sig (math calc "[M0]/1e17") 3) " x 10^17 N m, Mw = " (math round [Mw] 2)

Result

mentdb

isotropic = 0.2333, T N P axes = 3.505, 1.182, -4.687, double couple = 49.6 %, M0 = 4.22 x 10^17 N m, Mw = 5.68

Go further

Read the eigenvectors of matrix eigen to get the azimuth and plunge of the P and T axes.
Get strike, dip and rake of the two nodal planes from the T and P vectors.