Earthquake magnitudes

How big was it ?   That's what everyone wants to know about earthquakes.   Seismologists calculate earthquake magnitudes by measuring the amplitude of the wiggles recorded on seismometers and applying a correction factor to account for how far away the earthquake was.    The first person to do this was Charles Richter, he called the results of his calculations the Richter Magnitude.

In the UK seismologists use a formula that calculates the 'Local earthquake magnitude' (see https://academic.oup.com/gji/article/216/2/1145/5185119 )
 
     amp is seismogram displacement in nm (0.5x peak-peak) and  distance r is in  km

 ML=log(amp) + 1.11log(r) +0.00189r -1.16e^-0.2r -2.09

For the 8 August 2019 cornwall earthquake we can calculate the magnitude for each seismogram

station CCA1 magnitude Ml= 2.3
station RB30C magnitude Ml= 2.7
station RB5E8 magnitude Ml= 2.2
station RD93E magnitude Ml= 2.4
station R82BD magnitude Ml= 2.2
station R7FA5 magnitude Ml= 2.1

from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obspy import read, Stream
from obspy.geodetics import gps2dist_azimuth
from math import log10,exp

# station statistics
seislist=['CCA1','RB30C','RB5E8','RD93E','R82BD','R7FA5']
locations=['CM','FL','PZ','RD','RL','TS']
latitudes=[50.1867,50.1486,50.117,50.2344,50.2596,50.2609]
longitudes=[-5.2273,-5.0945,-5.535,-5.2384,-5.1027,-5.0434]

# Earthquakes' epicenter  BGS revised location
eq_lat=50.067
eq_lon=-5.2771

# set the data window
starttime=UTCDateTime("2019-08-08 16:52:05.10")
endtime=starttime+10

# start reading the data
# read in BGS seismometer at Carnmenellis
client=Client("IRIS")
waveform=client.get_waveforms('GB','CCA1','--','HHZ',starttime,endtime,attach_response= True)

# read in list of Raspberry Shakes by looping through the list from the second seismometer
client=Client(base_url='https://fdsnws.raspberryshakedata.com/')
for x in range (1,len(seislist)):
    waveform+=client.get_waveforms('AM',seislist[x],'00','EHZ',starttime,endtime,attach_response= True)

# work out the distance of each seismometer from the epicentre
for y in range(len(seislist)):
    waveform[y].stats["coordinates"] = {} # add the coordinates to your dictionary, needed for the section plot
    waveform[y].stats["coordinates"]["latitude"] = latitudes[y]
    waveform[y].stats["coordinates"]["longitude"] = longitudes[y]
    waveform[y].stats["distance"] = gps2dist_azimuth(waveform[y].stats.coordinates.latitude, waveform[y].stats.coordinates.longitude, eq_lat, eq_lon)[0]
    # Set the abbreviated name of the location in the network field for the plot title

    waveform[y].stats.network = locations[y]
    waveform[y].detrend(type='demean') ## remove any dc bias from traces

    # converting to velocity units using instrument sensitivity 
    # then and doing an integration to get displacement in metres
    waveform[y].remove_sensitivity()  #  data now in velocity units 
    waveform[y].integrate(method='spline')   #  data now in displacement units 
    waveform[y].filter('bandpass',freqmin=1.25,freqmax=18)  
    ## butterworth 4 pole bandpass filter using same parameters as BGS do for magnitude calculations

    # calculating the peak-peak amplitude of the trace
    amp_plus=max(waveform[y].data)
    amp_minus=min(waveform[y].data)
    waveform[y].stats["P_P_amplitude_nm"]=1000000000*(amp_plus-amp_minus)
    amp_nm=1000000000*(amp_plus-amp_minus)  ### for magnitude calculations amplitude in nanometres 
    r_km=0.001*waveform[y].stats["distance"]  ### distance must be in km 
    #  BGS magnitudes from https://academic.oup.com/gji/article/216/2/1145/5185119
    #  Extending local magnitude ML to short distances
    #  amplitude is displacement in nm (0.5x peak-peak) distance r is in  km
    #  ML=log(amp) + 1.11log(r) +0.00189r -1.16e^-0.2r -2.09
    #  note BGS use horizontal component data and measure the signal amplitude of the S/surface wave
    #  BGS also filter displacement traces with a 1.25-18Hz butterworth 4pole filter before calculating magnitudes

    waveform[y].stats["ML_magnitude"]=log10(amp_nm/2)+1.11*log10(r_km)+0.00189*r_km-1.16*exp(-0.2*r_km)-2.09
    print('station',seislist[y],'magnitude Ml=',round(waveform[y].stats["ML_magnitude"],1))

see  https://www.dentonseismo.co.uk  for more details of educational and citizen seismology and resources

Comments

Popular posts from this blog

lego seismometer

Barnsley Tyke quakes

Getting data from a BBC:microbit