Sunday, March 27, 2022

A Rough (and Silly) Calculation of the Energy Released by Earthquakes

A Facebook friend mentioned that he and his family were safe after the recent major earthquake in Taiwan, comparing it to the 1999 earthquake (which would be the Sept. 21 Chi-Chi earthquake, I guess).

I think it was his brother who jokingly asked whether 1991 was the year or the magnitude.

My tendency to go off on tangents was triggered, and I went looking for formulaic conversions of magnitude to energy released. 

(Yes, I know that the conversion is very much approximate. I also know that it's an exponential relationship, and that the theoretical limit for the earth is considered to be around 10. As I say, I tend to go off on tangents.)

I found some pages on earthquakes on the US Geological Survey site, including this page on magnitude ranges and their definitions

https://www.usgs.gov/programs/earthquake-hazards/earthquake-magnitude-energy-release-and-shaking-intensity
Partway down the page, I found this formula for relating moment magnitude to energy:

log E = 5.24 + 1.44M

That is the log (base ten) of the energy E is 1.44 times the moment magnitude plus 5.24.

Which means that the energy can be computed by the exponentiation 

105.24 + 1.44M

(You'll note that the units are not specified. I guess we can assume joules or ergs? More likely joules, 'though.)

So I got out my trusty Unix tool, bc, the multi-precision basic calculator, and proceeded to dig the definition of base ten exponentiation in terms of e out of my memory banks, somewhere.

The bc function

e(p)

is the the constant e raised to the power of the parameter p, the functional inverse of the natural logarithm. There is no predefined corollary for base ten, so we must define it ourselves, in terms of e.

10x = exln(10)

Invoke bc with the -l option on the command line:

bc -l

and we can define the power of ten function:

define power10(x) { return e(x*l(10)) }

But I knew I'd want a general exponentiation function, and I knew that bc loses precision, so I went to a little extra work:

define power(base,x) 
{ auto result;
  scale+=10;
  result = e(l(base)*x);
  scale-=10;
  return result
}

It's still a little rough, but close enough. Test it:

power(10,8)
99999999.999999999999999999999708633919

This shows one of the weaknesses of bc, but we can interpret that as very close to

100000000
Now we can define the earthquake energy function according to the definition above:

define energy(m) {return power(10,(5.24+1.44*m))}

And we can test it by running a loop to generate energies for reasonable magnitudes:

for ( i=0; i<=11; i+=0.5 ) { print i, ":\t", energy(i), "\n"; }

which produces the following:

0:    173780.082874937546699959956171209558
.5:    912010.839355909742120959407916159966
1.0:    4786300.923226383439220877393324318731
1.5:    25118864.315095801110850320677910503861
2.0:    131825673.855640710204737474230015115986
2.5:    691830970.918936487533684326974615305725
3.0:    3630780547.701013424673712123611850484194
3.5:    19054607179.632471826880141839831162456786
4.0:    99999999999.999999999999999999599371638788
4.5:    524807460249.772597364312157019920401025886
5.0:    2754228703338.166448631210659407094085487717
5.5:    14454397707459.275119314815458277888756899792
6.0:    75857757502918.376875376674665749727666287853
6.5:    398107170553497.250770252305085475893218954218
7.0:    2089296130854039.483122233735785788080932831027
7.5:    10964781961431850.131437136061343268093531268828
8.0:    57543993733715693.006203799395172142721740791878
8.5:    301995172040201619.863114517855660230944723734491
9.0:    1584893192461113485.202101373379733510028637417452
9.5:    8317637711026710061.666914027324396344905472398244
10.0:    43651583224016596746.383499610189244756061613362325
10.5:    229086765276777304572.408491985711096549196654929823
11.0:    1202264434617412905832.612715183424437635637611558644

It would be nice to get rid of the excess fractional digits, but bc is a little weak on formatting, so we won't work that hard.

But what can we check it against?

There is a table of magnitudes and energies on the same page. You can bring it into it's own browser window by clicking on it. Unfortunately, that table gives its results in terms of explosives. (Can we assume TNT?) Nor does it specify whether the magnitude shown is the moment magntude. And we are shown no easy way to convert. 

But using bc to get a ratio from the results of our formula to the table indicates a sort-of constant factor, so maybe we're on some sort of right track. 

energy(2)/56
2354029.89027929839651316918
energy(3)/1800
2017100.30427834079148539562
energy(4)/56000
1785714.28571428571428571428
energy(5)/1800000
1530127.05741009247146178369

Well, no, not really constant. We can hope, though, right?

A little thought gave me the idea that I could replicate the table with bc. 

Increasing the magnitude by two increases the energy by 103. In other words, squaring the base of the magnitude scale yields a thousand-fold increase in energy output. 

B2 <=>103

or

B <=> (103)1/2

so the base involved is related to the square root of 1000.

Continuing the sequence down to the missing magnitude 1 would give 1.8 kg of explosive. (And down one more to 0 would give 0.056 kg.)

logB(1.8) == 1

thus

1.8 == B1

That indicates a factor of 1.8.

Now we need a base 10 logarithm, or even a general logarithm, so we define that in bc:

define logg(base,val) { return l(val)/l(base)}

 And we can define the magnitude in terms of the energy:

define m(energy) { return logg(sqrt(1000),(energy/1.8)) }

A quick check shows we are close:

m(1800)
2.00000000000000000000
m(56)
.99527701460192956436
m(56000)
2.99527701460192956437

And we can define energy in terms of magnitude. Name it different from the above energy:

define enrgexp(m) { return power(sqrt(1000),m)*1.8 }

More quick checks:

enrgexp(2)
1799.999999999999999998982771182032
enrgexp(3)
56920.997883030827975931832400293898
enrgexp(4)
1799999.999999999999997965542364066557

Again, we're close. Maybe we wanted something more like 1.78 than 1.8 or something. But we'll call it close enough.

Let's generate the table:

for ( i = 0; i <= 11; i += 0.5 ) { print i, ":\t", enrgexp(i), "\n"; }
0:    1.800000000000000000000000000000
.5:    10.122143853426283447107688641472
1.0:    56.920997883030827975963999999956
1.5:    320.090293807006104220440146355010
2.0:    1799.999999999999999998982771182032
2.5:    10122.143853426283447101968343457696
3.0:    56920.997883030827975931832400293898
3.5:    320090.293807006104220259254648782792
4.0:    1799999.999999999999997965542364066557
4.5:    10122143.853426283447096248045442524754
5.0:    56920997.883030827975899664800630542382
5.5:    320090293.807006104220078362942553947890
6.0:    1799999999.999999999996948313546099836862
6.5:    10122143853.426283447090527747427352572037
7.0:    56920997883.030827975867497200967185895594
7.5:    320090293807.006104219897471236325102741328
8.0:    1799999999999.999999995931084728133115816211
8.5:    10122143853426.283447084807449412180389148809
9.0:    56920997883030.827975835329601303829406636110
9.5:    320090293807006.104219716579530096257591698311
10.0:    1799999999999999.999994913855910166394770266991
10.5:    10122143853426283.447079087151397008206259475075
11.0:    56920997883030827.975803162001640472917677075752

This is very close to the table. Let's compare the two calculations:

0:    1.800000000000000000000000000000 --    96544.49048607641483331108
.5:    10.122143853426283447107688641472 --    90100.56096438499849155669
1.0:    56.920997883030827975963999999956 --    84086.73602423378680358317
1.5:    320.090293807006104220440146355010 --    78474.3080345974612309003\
6
2.0:    1799.999999999999999998982771182032 --    73236.485475355950113784\
42
2.5:    10122.143853426283447101968343457696 --    68348.26504513230745791\
163
3.0:    56920.997883030827975931832400293898 --    63786.31230538237508595\
543
3.5:    320090.293807006104220259254648782792 --    59528.8502909781343005\
2010
4.0:    1799999.999999999999997965542364066557 --    55555.555555555555555\
61834
4.5:    10122143.853426283447096248045442524754 --    51847.46115538839130\
740462
5.0:    56920997.883030827975899664800630542382 --    48386.86610867114652\
844766
5.5:    320090293.807006104220078362942553947890 --    45157.2508980055136\
9680788
6.0:    1799999999.999999999996948313546099836862 --    42143.198612732431\
59750293
6.5:    10122143853.426283447090527747427352572037 --    39330.32135467432\
837797306
7.0:    56920997883.030827975867497200967185895594 --    36705.19155597755\
591727517
7.5:    320090293807.006104219897471236325102741328 --    34255.2778811934\
5368453639
8.0:    1799999999999.999999995931084728133115816211 --    31968.885407619\
82944796326
8.5:    10122143853426.283447084807449412180389148809 --    29835.09979834\
737393915950
9.0:    56920997883030.827975835329601303829406636110 --    27843.73520151\
512694460159
9.5:    320090293807006.104219716579530096257591698311 --    25985.2856270\
6656631018588
10.0:    1799999999999999.999994913855910166394770266991 --    24250.87956\
889810930361491
10.5:    10122143853426283.447079087151397008206259475075 --    22632.2376\
5578404058365912
11.0:    56920997883030827.975803162001640472917677075752 --    21121.6331\
2892006645282710

Ick. Definitely not constant. We need more research. 

Oh, well, let's see what these two give us for magnitude 1999. (Expect the computer to be thinking for several seconds, or, if it's an older computer, for several minutes.) 

According to the table:

enrgexp(1999)
56920997883030827943828567936264122901386094362240781432147822485612\
46108160721544492433023459350568540095430581546108778254387446161457\
20070148007795573394566997681517496549855105967989636358346558093359\
47073049417418282438320742899890844575386635960109440927286954346859\
26664139131795793863000643693800633418743164057846984015963204398775\
60449145286650344364978135609191881503644766064649294148861402387152\
39668774258714917538971629308850707288043396599541120218142796686010\
43452996292621622704149840627054832081858481088646690119679417531046\
56298402816476046016168622753716951661114147884630225967967062033982\
25898508604818718666594848139667173044074122438610000765611500984899\
77239324279940490269201971705103863769451211341837543162365120094953\
96372566379175635507597956561845806880071329453040286902000420349163\
85556066776612019516479819200708279918388326184232205610167337240035\
94968674438164951105098989018168645726827415188515197163180391272780\
96854853137878879070351891132065229334557014601509853950423709113692\
71504168355398757242593228019984379655528838024245708723589897128072\
06498876979242507707136599994074216464282769916020834708110157811783\
92097264654920969115336939219106453572424091851472077880805266451310\
78892060527925402753269382486372395414335463983902503563715677564875\
20857056570243196764478479653237433742689869602654557728934702873455\
41058528214172546741565064845323537805309209583322968601762463253981\
65768834295654907389088071030176830969776690136997679582726069548353\
28773902457115082080246159058128134088413409687693428856748628952721\
62561036856811787872061563354575980408785346723568831127928828263324\
36833042000703289399191833063516040190581018961872031891034612931273\
32627551884393379281377030408790506005097096839816761861055853106847\
40273913392304829522168642644514093315208401226560551081020631661869\
44669362924382340029541873299612457575267645299410112616845610609712\
57155632112179595072545779921282292174073071106414165441012581922077\
58240912573195894271374086191892707980071349386026984672449337392590\
71811036521144830194808428047529184528351324268579670472146830818718\
03658914353990832099189312973584391403704614271720373891068337353002\
19760291339060945982644854953925767185052858435916806267544829757190\
35023347003009930833854204312056981331410705410406981823893847010311\
65966359687128923907759705435892521150916017836150226473121465388121\
65244041088025014105721061272942790251956828459146625113216736414629\
01255285568149140085256112880790465496296451782787238210648925532468\
91228244517803214061916621548124626660136148606174974948739406508499\
75958850823187858133382225870765447783329396941523348286236042738989\
97271611205467743793586953807091177880964118682802458622738494994943\
66638238646583547730809130783602166998836765721270540463411437221582\
20761071981983598060411854943722750584070642863371824436238575540160\
32012236310452853863877998393972886045626427628828990733337222861559\
33076566577984941602988965669084812571961771795594798948084265337304\
7091990.972722246071529279958378058522

According to the formula:

energy(1999)
63095734448019324943436013595952198806106263004472978130750283914739\
85646172253254023959427080334085792471330807169010078706358182486189\
60033689864917960134967706939099353312486237647979666174700395609250\
42451159872490899438738395885908500940174729834832255176815048574805\
61661450752614317546358515330045104383717404214399019027819009914096\
64357009583667771854640635722153443357824965277318111125501776090051\
96062103303822867056576883905202092068873898904542385677847712405612\
77565656173618332709652305654968664984233313775929647720898582312099\
70156272622210587374188362416663183736841266890733504386411018565642\
08304922038813848527770051224186152002307571087935067658324537519247\
65407671367387523766404242701715874108621296255965529985298096778568\
58248551995469123897978822819952278506945693873982434642215771322665\
69652926166579216908501746332112830698445697667370000842277897407994\
03398930095661631038032448414964667985529179433521450906452665656153\
50791195318492017696436071528350459093931459644638582402920662216460\
76747836124122291219107303624228131486076883477323983932150580507123\
00457162434527509114584761617503767004925432456579888693189401934784\
38482698018809664859305044757007774590548167569512868829690186401774\
27532653628927265722925135823326887893322744136616091407587392180803\
65179668096605267790786219685059270845270549609370914640375194156385\
22209220950613357145003579842368213095807264933735537752554310524517\
28580202244636638758277110505875140283236900456123668156484687844417\
38553837772926253755441635181187918998242081233055931794909345873601\
79282653801644611750682054117429315994400974514247776958824343456877\
53226461800417288452729321436837612317061388523209381562589480514804\
89799607307813720614394034642302784994662135916160246580689941278950\
52686494945485061168988117965157865756367116685076682781272561899033\
49195531987729238602724220712896587312066519658647207530726171012380\
77943681712018379283403905515367085045196938214489856330547269669281\
67928782077756768955437584316044790771414158399097879703387896284391\
99483961037533641447959247048894834920911650159157633405470028053377\
61781655053832962677338524122974477967065186981728961367238358136305\
17523483808265012311416078069441513670245974603206188308861526642599\
47597753125583196306804658491680492337057681325060525106325330483736\
88462800730318900305747536613840534202550304481866710836074367345455\
99288067975645066852701427952293314445755720927898179618294864028322\
65706420727926722366060858167193117503791685342349860784112828071626\
30092883969265869257068422309693548490502597994378924632947416370609\
37180882800242053789139556220177794975223407421317091737933310766707\
30032550700165344079288743950437998641164201591698807189441678574528\
59496472666250949210911796728354955093122911558869350963905286175742\
06211260186702048098154989247956439941034888415126044388808567869251\
6943222441701152025362619807.423572390462876603381595102367

Having a hard time even counting the digits on those? Copy them into a text edit document and delecte the backslash-newlines and you can use the character count tool or even just the column count indicator of your text editor to count that. (Or you can use your shell and the tools tr and wc.) 

Remember to delete the fractional digits at the ends.

enrgexp(1999) 56920997883030827943828567936264122901386094362240781432147822485612461081607215444924330234593505685400954305815461087782543874461614572007014800779557339456699768151749654985510596798963635834655809335947073049417418282438320742899890844575386635960109440927286954346859266641391317957938630006436938006334187431640578469840159632043987756044914528665034436497813560919188150364476606464929414886140238715239668774258714917538971629308850707288043396599541120218142796686010434529962926216227041498406270548320818584810886466901196794175310465629840281647604601616862275371695166111414788463022596796706203398225898508604818718666594848139667173044074122438610000765611500984899772393242799404902692019717051038637694512113418375431623651200949539637256637917563550759795656184580688007132945304028690200042034916385556066776612019516479819200708279918388326184232205610167337240035949686744381649511050989890181686457268274151885151971631803912727809685485313787887907035189113206522933455701460150985395042370911369271504168355398757242593228019984379655528838024245708723589897128072064988769792425077071365999940742164642827699160208347081101578117839209726465492096911533693921910645357242409185147207788080526645131078892060527925402753269382486372395414335463983902503563715677564875208570565702431967644784796532374337426898696026545577289347028734554105852821417254674156506484532353780530920958332296860176246325398165768834295654907389088071030176830969776690136997679582726069548353287739024571150820802461590581281340884134096876934288567486289527216256103685681178787206156335457598040878534672356883112792882826332436833042000703289399191833063516040190581018961872031891034612931273326275518843933792813770304087905060050970968398167618610558531068474027391339230482952216864264451409331520840122656055108102063166186944669362924382340029541873299612457575267645299410112616845610609712571556321121795950725457799212822921740730711064141654410125819220775824091257319589427137408619189270798007134938602698467244933739259071811036521144830194808428047529184528351324268579670472146830818718036589143539908320991893129735843914037046142717203738910683373530021976029133906094598264485495392576718505285843591680626754482975719035023347003009930833854204312056981331410705410406981823893847010311659663596871289239077597054358925211509160178361502264731214653881216524404108802501410572106127294279025195682845914662511321673641462901255285568149140085256112880790465496296451782787238210648925532468912282445178032140619166215481246266601361486061749749487394065084997595885082318785813338222587076544778332939694152334828623604273898997271611205467743793586953807091177880964118682802458622738494994943666382386465835477308091307836021669988367657212705404634114372215822076107198198359806041185494372275058407064286337182443623857554016032012236310452853863877998393972886045626427628828990733337222861559330765665779849416029889656690848125719617717955947989480842653373047091990.972722246071529279958378058522

That's 3014 digits to the left of the decimal point. Order of 103014.

energy(1999)
6309573444801932494343601359595219880610626300447297813075028391473985646172253254023959427080334085792471330807169010078706358182486189600336898649179601349677069390993533124862376479796661747003956092504245115987249089943873839588590850094017472983483225517681504857480561661450752614317546358515330045104383717404214399019027819009914096643570095836677718546406357221534433578249652773181111255017760900519606210330382286705657688390520209206887389890454238567784771240561277565656173618332709652305654968664984233313775929647720898582312099701562726222105873741883624166631837368412668907335043864110185656420830492203881384852777005122418615200230757108793506765832453751924765407671367387523766404242701715874108621296255965529985298096778568582485519954691238979788228199522785069456938739824346422157713226656965292616657921690850174633211283069844569766737000084227789740799403398930095661631038032448414964667985529179433521450906452665656153507911953184920176964360715283504590939314596446385824029206622164607674783612412229121910730362422813148607688347732398393215058050712300457162434527509114584761617503767004925432456579888693189401934784384826980188096648593050447570077745905481675695128688296901864017742753265362892726572292513582332688789332274413661609140758739218080365179668096605267790786219685059270845270549609370914640375194156385222092209506133571450035798423682130958072649337355377525543105245172858020224463663875827711050587514028323690045612366815648468784441738553837772926253755441635181187918998242081233055931794909345873601792826538016446117506820541174293159944009745142477769588243434568775322646180041728845272932143683761231706138852320938156258948051480489799607307813720614394034642302784994662135916160246580689941278950526864949454850611689881179651578657563671166850766827812725618990334919553198772923860272422071289658731206651965864720753072617101238077943681712018379283403905515367085045196938214489856330547269669281679287820777567689554375843160447907714141583990978797033878962843919948396103753364144795924704889483492091165015915763340547002805337761781655053832962677338524122974477967065186981728961367238358136305175234838082650123114160780694415136702459746032061883088615266425994759775312558319630680465849168049233705768132506052510632533048373688462800730318900305747536613840534202550304481866710836074367345455992880679756450668527014279522933144457557209278981796182948640283226570642072792672236606085816719311750379168534234986078411282807162630092883969265869257068422309693548490502597994378924632947416370609371808828002420537891395562201777949752234074213170917379333107667073003255070016534407928874395043799864116420159169880718944167857452859496472666250949210911796728354955093122911558869350963905286175742062112601867020480981549892479564399410348884151260443888085678692516943222441701152025362619807.423572390462876603381595102367

And that one is 2885 digits to the left of the decimal point. Order of 102885.

What does that even mean?

Checking the energy output of the sun, this is roughly in the range of the total energy output of the sun for several seconds -- concentrated on little old earth.

Heh. Quite a silly exercise. Tools like bc can be a little dangerous, in terms of wasting time. But it does give one an appreciation of how meaningless magnitude numbers become after a certain point.