Seismology and Isostasy#

IN THE SPACE BELOW, WRITE OUT IN FULL AND THEN SIGN THE HONOR PLEDGE:

“I pledge my honor that I have not violated the honor code during this examination.”

PRINT NAME:

If a fellow student has contributed significantly to this work, please acknowledge them here:

Peer(s):

Contribution:

By uploading this assignment through Canvas, I sign off on the document below electronically.

Downloads from the Internet#

We will need some files from the internet that will be used in this assignment. If you get errors in the following cell block, it probably means the machine (or the node in a cluster) you are calling this command from does not have direct internet access.

Note: If you are on the Adroit cluster in Princeton, please make sure you are within a Jupyter on Adroit Vis session since this node has internet access. You can use any other session for the remaining cells in this notebook; Jupyter for classes is suggested since a dedicated class environment (GEO203) with all necessary modules has been created for this course.

import os
import requests
from pathlib import Path

# Define the list of remote files to retrieve\n",
remote_urls = ['https://portal.globalseismology.org/files/DATA/GLOBALTOPO/ETOPO1_Ice_g_gmt4.grd']
path = '../PS_Plate_Tectonics/Files/'

# loop through all remote URL and download the files
for remote_url in remote_urls:
    # Define the local filename to save data
    # get the trailing filename after the last / of the path
    local_file = path + os.path.basename(remote_url)
    
    if not Path(local_file).is_file():
        # mke the directory if it does not exist
        if not os.path.isdir(path): os.makedirs(path)
    
        # Download remote and save locally
        r = requests.get(remote_url, allow_redirects=True)
        open(local_file, 'wb').write(r.content)

Part II: Seismic stations, seismographs, and seismograms#

Seismographs are instruments at seismic stations that measure ground motion. The data is displayed as seismograms. This data is of fundamental importance to the geosciences. By studying earthquake-generated ground motion we can learn about 1) the nature of the earthquake itself (the source of the motion) and 2) the material the earthquake waves travel through (Earth’s interior). This is because the observed vibration (in the time domain t) is the result of the effects of earthquake source, interior structure and how the instrument senses the vibration. These effects can be most easily modeled by considering each frequency (\(\omega\)) of the waveform in turn, and using their products (i.e. frequency-domain convolution per mathematical definition) to derive the waveform.

Figure: Relationship between measurement and contribution from instrument, structure and source. The United States Geologic Survey estimates that several million earthquakes occur each year. Many go undetected because they hit remote areas without good seismic station coverage, or have very small magnitudes. But various seismic networks locate something like 20,000 earthquakes a year. We will focus on two large earthquakes since they caused the whole Earth to ring like a bell, causing waveforms to be recorded at various stations well above the noise level of the instrument. Tohoku earthquake USGS Link The March 11, 2011, M 9.1 Tohoku earthquake, which occurred near the northeast coast of Honshu, Japan, on the subduction zone plate boundary between the Pacific and North America plates. The location, depth (about 25 km), and focal mechanism solutions of the March 11th earthquake are consistent with the event having occurred on the subduction zone plate boundary. Modeling of the rupture of this earthquake indicates that the fault moved as much as 50–60 m, and slipped over an area approximately 400 km long (along strike) by 150 km wide (in the down-dip direction). The Japan Trench subduction zone has hosted nine events of \(m_w\) 7+ since 1973.

Mini-tutorial#

Reading the CMT solution file for an earthquake

The focal mechanisms that are usually given to convey this configuration are a way of plotting the three dimensional fault information in two dimensions (they’re lower hemisphere stereographic projections). The figure below is an attempt to show how a simple block model of the first motions of 5 common types of earthquake look when projected onto a lower hemisphere, and then when the pattern on that focal sphere is projected onto the horizontal to produce a 2D focal mechanism.

Figure: Types of faulting and their expression as focal mechanism of the earthquake. The Global Centroid-Moment-Tensor (CMT) Project has measured the moment tensors of all earthquakes with magntiudes \(m_w\) 5.5+ since 1976. The CMT project has been continuously funded by the National Science Foundation since its inception. Let us read the CMT solution of the Tohoku earthquake in the CMTSOLUTION format, available from their webpage.

import numpy as np
import os
import subprocess
from multiprocessing import Pool
import glob
import errno
import cmath
from math import atan2
import pandas as pd
import pdb
from obspy import UTCDateTime,read_events,read,read_inventory
from obspy.geodetics import gps2dist_azimuth

#########################################################

def readcmtsolution(file='CMTSOLUTION'):
    if not os.path.isfile(file): raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), file)
    cmt={}
    f = open(file,'r')
    line=f.readline().split()
    cmt['source']=line[0];cmt['iyr']=int(line[1]);cmt['imo']=int(line[2]);cmt['ida']=int(line[3])
    cmt['iho']=int(line[4]);cmt['imi']=int(line[5]);cmt['fsec']=float(line[6]);cmt['eplat']=float(line[7]
)
    cmt['eplong']=float(line[8]);cmt['depth']=float(line[9]);cmt['mb']=float(line[10]);cmt['ms']=float(line[11])
    cmt['region']=' '.join(line[12:])

    cmt['ityphdur']=2
    cmt['itypso']=0
    cmt['time'] = UTCDateTime(cmt['iyr'],cmt['imo'],cmt['ida'],cmt['iho'],cmt['imi'],cmt['fsec'])

    # read cmts
    xm = np.zeros(6)
    while line:
        line = f.readline()
        if line[0:10] == 'event name':
            cmt['name']=line[12:].split()[0]
        elif line[0:10] == 'time shift':
            cmt['tcmt']=float(line[12:].split()[0])
        elif line[0:13] == 'half duration':
            cmt['hdur']=float(line[15:].split()[0])
        elif line[0:8] == 'latitude':
            cmt['epla']=float(line[10:].split()[0])
        elif line[0:9] == 'longitude':
            cmt['eplo']=float(line[11:].split()[0])
        elif line[0:5] == 'depth':
            cmt['depth']=float(line[7:].split()[0])
        #Double  couple
        elif line[0:3] == 'Mrr':
            xm[0]=float(line[5:].split()[0])
        elif line[0:3] == 'Mtt':
            xm[1]=float(line[5:].split()[0])
        elif line[0:3] == 'Mpp':
            xm[2]=float(line[5:].split()[0])
        elif line[0:3] == 'Mrt':
            xm[3]=float(line[5:].split()[0])
        elif line[0:3] == 'Mrp':
            xm[4]=float(line[5:].split()[0])
        elif line[0:3] == 'Mtp':
            xm[5]=float(line[5:].split()[0])
        # Single force
        elif line[0:3] == 'CFr':
            cmt['itypso']=1
            xm[1]=float(line[5:].split()[0])
        elif line[0:3] == 'CFt':
            cmt['itypso']=1
            xm[2]=float(line[5:].split()[0])
        elif line[0:3] == 'CFp':
            cmt['itypso']=1
            xm[3]=float(line[5:].split()[0])
        # Single force
        elif line[0:3] == 'SFr':
            cmt['itypso']=2
            xm[1]=float(line[5:].split()[0])
        elif line[0:3] == 'SFt':
            cmt['itypso']=2
            xm[2]=float(line[5:].split()[0])
        elif line[0:3] == 'SFp':
            cmt['itypso']=2
            xm[3]=float(line[5:].split()[0])

    if cmt['itypso'] > 1:
        xm[3]=0.
        xm[4]=0.
        xm[5]=0.
    cmt['xm']=xm
    
    # Moment tensor
    M = np.array([[xm[0], xm[3], xm[4]], [xm[3], xm[1], xm[5]], [xm[4], xm[5], xm[2]]])
    
    # Scalar moment
    M0 = np.sqrt(np.sum(M ** 2) / 2)
    cmt['m0'] = M0
    
    # Moment magnitude
    Mw = np.round(np.log10(M0) / 1.5 - 10.73, 1)
    cmt['mw'] = Mw

    return cmt

# read the CMT solution file of Tohoku earthquake
cmt = readcmtsolution('Files/CMTSOLUTION_TOHOKU')
print(cmt)

{'source': 'PDEW', 'iyr': 2011, 'imo': 3, 'ida': 11, 'iho': 5, 'imi': 46, 'fsec': 24.68, 'eplat': 38.32, 'eplong': 142.351, 'depth': 20.0, 'mb': 7.2, 'ms': 9.0, 'region': 'TOHOKU', 'ityphdur': 2, 'itypso': 0, 'time': UTCDateTime(2011, 3, 11, 5, 46, 24, 680000), 'name': '201103110546A', 'tcmt': 69.76, 'hdur': 70.0, 'epla': 37.52, 'eplo': 143.05, 'xm': array([  1.73000000e+29,  -2.81000000e+28,  -1.45000000e+29,
         2.12000000e+29,   4.55000000e+29,  -6.57000000e+28]), 'm0': 5.3118480305822004e+29, 'mw': 9.0999999999999996}

from obspy.imaging.beachball import beachball
fig = beachball(cmt['xm'])

TO DO#

Question 3.1 Based on the beachball focal mechanism above, what type of faulting caused the Tohoku earthquake? Was the fault shallowly dipping at a low angle from the surface (dip < 30\(^{\circ}\)) or steeply dipping at a high angle (dip > 60\(^{\circ}\))? Which plate was the overriding plate? Answer:

Section A: Signal processing of vibrations at a single station#

Modern seismometers measure ground motion in three mutually perpendicular directions:

North-South (LHN or BHN component). This records the north-south horizontal component of vibration. “Up” on this graph means the ground moved north. “Down” on this graph means the ground moved south.
East-West (LHE or BHE component). This records the east-west horizontal component of vibration. “Up” on this graph means the ground moved east. “Down” on this graph means the ground moved west.
Vertical (LHZ or BHZ component). This records the vertical component of vibration. “Up” on this graph the means the ground moved up. “Down” on this graph means the ground moved down.

The prefix (L-long period, or B-broadband) indicates the sampling rate of seismic data stream. When we wish to analyse signals that are finely spaced, we use the broadband data steam (i.e. think of Nyquist frequency). Sometimes, the horizontal orientations of the sensor in the instrument are not perfectly aligned with the N-S and E-W directions, even if the vertical component is perfectly aligned. In these cases, the horizontal components are denoted as LH1 and LH2 (or BH1 and BH2). We typically analyze the motion along the propagation direction of the wavefront, which gives us clearer information on individual seismic phases (as you will see below). Note that finding the propagation direction from earthquake to a station requires the coordinates and depth of the hypocenter, as denoted by the CMT focal mechanism above.

Figure: Rotation of vibrations recorded on a seismogram to radial, transverse and vertical components along the propagation direction. Let us perform this calculation for the vibrations recorded on the long-period components at Ascension Island (station code: ASCN) in the South Atlantic Ocean from the Tohoku earthquake.

st = read('Files/waveforms/II.ASCN.*.LH*.mseed') # Ascension Island, distance=144
print(st)

3 Trace(s) in Stream:
II.ASCN.00.LH1 | 2011-03-11T05:46:25.069500Z - 2011-03-11T15:46:24.069500Z | 1.0 Hz, 36000 samples
II.ASCN.00.LH2 | 2011-03-11T05:46:25.069500Z - 2011-03-11T15:46:24.069500Z | 1.0 Hz, 36000 samples
II.ASCN.00.LHZ | 2011-03-11T05:46:25.069500Z - 2011-03-11T15:46:24.069500Z | 1.0 Hz, 36000 samples

st.plot();

We want to convert from counts (e.g. voltage) reported by the seismometer above to ground dispacement (or velocity or accelaration), which is the value we are interested in. This is done by considering the instrument response of the instrument. For this station (ASCN), the response belows shows that all vibrations at periods roughly 1-3000 seconds will be recorded with enough amplitude counts.

stationdir = 'Files/stations/'
inv = read_inventory('%s%s.%s.xml' % (stationdir, st[0].stats.network, st[0].stats.station))
inv.plot_response(0.001, output="DISP");

We also wish to focus on vibrations in the period range 20-200 seconds to start with. This done by band-pass filtering the seismogram below:

periods = [20,200]

stat_lat = inv.select(station=st[0].stats.station)[0][0].latitude
stat_lon = inv.select(station=st[0].stats.station)[0][0].longitude
evt_lat = cmt['epla']
evt_lon = cmt['eplo']

# compute back azimuth (check sign, later!):
dist, back_az, az = gps2dist_azimuth(evt_lat, evt_lon, stat_lat, stat_lon)

# convert distance to degrees
dist = (dist / 6371000) / np.pi * 180

# select the traces from a station
st3 = st.copy()

# remove instrument response, filter, and rotate
pre_filt = [0.8/periods[1] ,1/periods[1], 1/periods[0], 1.25/periods[0]]
st3 = st3.remove_response(inventory=inv, output="DISP",pre_filt=pre_filt).taper(0.05)
if not st3.select(channel='LHE'): st3 = st3.rotate(method="->ZNE", inventory=inv,back_azimuth=back_az)
st3 = st3.rotate(method="NE->RT", inventory=inv,back_azimuth=back_az)
st3.plot();

import holoviews as hv
from holoviews.operation.datashader import datashade
from holoviews.streams import PlotSize
import datashader as ds
from holoviews import opts
PlotSize.scale=2
hv.extension('plotly')

import plotly
import plotly.graph_objects as go
import plotly.express as px
plotly.offline.init_notebook_mode() 

from obspy.taup import TauPyModel
from obspy.geodetics.base import gps2dist_azimuth
import warnings
warnings.filterwarnings("ignore")

Section B: Horizontal versus Vertical Motions#

Task I: At a single station#

duration_h = 1.5
st3 = st3.slice(endtime=st3[0].stats.starttime + 3600*duration_h)

# grab the time series data for processing
scalefactor = 0.5
times = np.array([])
signals = np.array([])
component_indicator = []
components = ["T", "Z"]
for ii, component in enumerate(components):
    trace = st3.select(component=component)[0]
    t = trace.times()
    x = trace.data
    x_scale = np.max(np.abs(x)) * scalefactor
    times = np.hstack((times, t/60))       # convert to minutes
    signals = np.hstack((signals, 2*ii + x / x_scale))
    component_indicator += [component] * t.size
        
# convert to dataframe
data = dict()
cols = []

data['time (minutes)'] = times 
data['signal (normalized)'] = signals
data['component'] = component_indicator
data['station'] = ['%s-%s' % (st3[0].stats.network, st3[0].stats.station)] * times.size
data['distance'] = [dist] * times.size
data['azimuth'] = [az] * times.size
dfm = pd.DataFrame(data)

fig_title = 'M%.1f %s lat=%.2f lon=%.2f dep=%.2f km band=%d-%d s station=%s-%s' % (cmt['mw'], 
                                                                                   cmt['region'].capitalize(), 
                                                                                   cmt['epla'], 
                                                                                   cmt['eplo'], 
                                                                                   cmt['depth'],
                                                                                   periods[0], 
                                                                                   periods[1],
                                                                                   st3[0].stats.network,
                                                                                   st3[0].stats.station)
fig_1 = px.line(dfm, x='time (minutes)', y='signal (normalized)', color='component', 
                hover_data={'station': ':.s', 'distance': ':.2f', 'azimuth':':.2f'}, 
                title=fig_title)
fig = go.Figure(data=fig_1,layout=fig_1.layout)
plotly.offline.iplot(fig)

TO DO#

Question 3.2 Identify the first seismic wave arrival at your station on the Z component. This is the first P-wave. Note that you should consider when the seismograph first starts to respond to the P-wave – i.e., the “shoulder” of the wave, or where trace starts to climb or fall towards the first peak or trough. Report your P-wave travel time pick in minutes and seconds.

Answer:

Question 3.3 Identify the first S-wave arrival at your station. Hint: how would the relative horizontal vs. vertical motion differ for an S wave compared to a P wave? There may also be a change in wavelength and/or amplitude; S-waves look ‘different” than P-waves. But there is likely more uncertainty about when to “pick” the arrival time of the first S-wave. This is because by this time, the station is still responding to the passage of P and other waves; and the S-waves are superimposed on these other waves. This is part of the “art” of seismology – and also partly why we use data from many stations to analyze seismic waves. Report your S-wave travel time pick in minutes and seconds.

Answer:

Question 3.4 Identify the surface waves on the seismogram. These are the waves with the strongest vibrations (largest amplitudes). Provide your travel time picks for the Rayleigh and Love wave packets on Z and T components, respectively. Note: Picking a travel time of a surface wave is complicated and depends strongly on the frequency of the surface wave, due to a property called dispersion. For this task, simply report the time of the first shoulder from the left (earliest in time). Which of the two waves (Love or Rayleigh) is arriving sooner at the station?

Answer:

Task II: Record Section of the USArray Transportable Array#

A single seismometer can give us some information about the earthquake and interior structure. In fact, this is the idea behind deploying a seismometer on the Moon and Mars. However, arrays of seismometer can provide unprecedented insights into geological processes.

USArray is a 15-year program to place a dense network of permanent and portable seismographs across the continental United States. The seismographs record local, regional, and distant (teleseismic) earthquakes. This mega-project is funded by the US National Science Foundation and we will use data from a part of this project called the Transportable Array.

Figure: Installation plan of the transportable array of USArray.

USArray has a large traveling network of 400 high-quality, portable seismographs that are being placed in temporary sites across the United States in a rolling fashion. Station spacing is ~70 km (42 mi). Transportable Array data are extremely useful for mapping the structure of Earth’s uppermost 70 km. The array was initially deployed in the westernmost United States. Unless adopted and made into a permanent installation, after 18-24 months, each instrument is picked up and moved to the next carefully spaced array location to the east. When completed, over 2000 locations will have been occupied during this program.

Note that since we are looking at the data from Tohoku earthquake in 2011, we are using data from the Transportable Array as it was installed at that time. This means the stations that we are analyzing lie on a grid from Michigan to Louisiana.

from scipy import interpolate

def plottraces(cmt, 
               st, 
               stationdir, 
               periods=[20,100],
               component="T", 
               phaselist=["P", "S", "PcP", "ScS"],
               plottype="holoviews",
               duration_h=1.0,
               scalefactor=1.0):
    '''
    Make holoview objects for plotting traces against epicentral distance
    
    Input:
    -----
    cmt:            CMT Solution
    st:             Obspy stream object containing 3 components seismograms
    stationdir:     directory to station XML files for the streams
    periods:        [min, max] periods for bandpass filter
    phaselist:      list of phases to be plotted
    plottype:       either "holoviews" or "datashader"
    duration:       the number of hours of seismograms to be plotted
    scalefactor:    factor to reduce the amplitude of the seismogram
    
    Output:
    ------
    overlay_trace:  holoview object for trace plots
    overlay_phase:  holoview object for phase plots
    options:        options used by datashader
    
    These outputs will be used in the next command to create a plot:
    ---------------------------------------------------------------
    datashade(ndoverlay * overlay_phase, cmap=["black"], cnorm='linear', 
              aggregator=ds.count(), line_width=8).opts(opts)
    
    '''
        
    ########################################################
    # READ, PREOCESS, AND GATHER TRACES
    ########################################################
    signals = []
    dists = []

    # set the time for all trace to be the same
    t_trace = st[0].slice(endtime=st[0].stats.starttime + 3600*duration_h).times()

    for tr in st.select(channel='LHZ'):
        inv = read_inventory('%s%s.%s.xml' % (stationdir, tr.stats.network, tr.stats.station))
        stat_lat = inv.select(station=tr.stats.station)[0][0].latitude
        stat_lon = inv.select(station=tr.stats.station)[0][0].longitude
        evt_lat = cmt['epla']
        evt_lon = cmt['eplo']

        # compute back azimuth (check sign, later!):
        dist, back_az, _ = gps2dist_azimuth(evt_lat, evt_lon, stat_lat, stat_lon)

        # convert distance to degrees
        dist = (dist / 6371000) / np.pi * 180
        dists += [dist]

        # select the traces from a station
        st3 = st.select(station=tr.stats.station).copy()

        # remove instrument response, filter, and rotate
        pre_filt = [0.8/periods[1] ,1/periods[1], 1/periods[0], 1.25/periods[0]]
        st3 = st3.slice(endtime=st3[0].stats.starttime + 3600*duration_h).remove_response(inventory=inv, 
                            output="DISP",pre_filt=pre_filt).taper(0.05)
        if not st3.select(channel='LHE'):
            st3 = st3.rotate(method="->ZNE", inventory=inv,back_azimuth=back_az)
        st3 = st3.rotate(method="NE->RT", inventory=inv,back_azimuth=back_az)

        # grab the time series data for processing
        trace = st3.select(component=component)[0]
        t = trace.times()
        x = trace.data

        # interpolate to the same time
        f = interpolate.interp1d(t, x, kind='linear', fill_value='extrapolate')
        x = f(t_trace)
        x_scale = np.max(np.abs(x)) * scalefactor
        signals += [dist + x / x_scale]
        
        
    # convert to dataframe
    data = dict()
    cols = []

    for ii in range(len(signals)):
        data["trace_%03d" % ii] = signals[ii]
        cols += ["trace_%03d" % ii]
    data['time'] = t / 60                                 # convert to minutes
    dfm = pd.DataFrame(data)
        
    ########################################################
    # CALCULATE TRAVEL TIMES
    ########################################################
    dist_min = np.min(dists)
    dist_max = np.max(dists)

    distances = np.arange(np.floor(dist_min), np.ceil(dist_max) + 1, 0.5)
    traveltime_phases = np.ones((np.size(distances), len(phaselist))) * np.nan

    model = TauPyModel(model="ak135")
    for ii, dist in enumerate(distances):
        for jj, phase in enumerate(phaselist):
            arrivals = model.get_travel_times(source_depth_in_km=cmt['depth'],
                                              distance_in_degree=dist,
                                              phase_list=[phase])
            if arrivals:
                traveltime_phases[ii, jj] = arrivals[0].time
                
    # convert to dataframe
    traveltime = dict()
    for jj, phase in enumerate(phaselist):
        traveltime[phase] = traveltime_phases[:,jj] / 60  # convert to minutes
    traveltime['distance'] = distances
    dft = pd.DataFrame(traveltime)
    
    ########################################################
    # PLOT
    ########################################################
    options = hv.opts.RGB(width=800, height=800, xlabel='Time (minutes)', 
                       title='M%.1f %s lat=%.2f lon=%.2f dep=%.2f km comp=%s band=%d-%d s, phase=%s' % (cmt['mw'], 
                                                                                                        cmt['region'].capitalize(), 
                                                                                                        cmt['epla'], 
                                                                                                        cmt['eplo'], 
                                                                                                        cmt['depth'], 
                                                                                                        component, 
                                                                                                        periods[0], 
                                                                                                        periods[1],
                                                                                                        phaselist))
    
    
    if plottype == "datashader":
        overlay_trace = hv.NdOverlay({c:hv.Curve((dfm['time'], dfm[c]), vdims = ['Distance (degree)'], label='trace') for c in cols})
        if len(phaselist) > 0: 
            overlay_phase = hv.NdOverlay({phase:hv.Curve((dft[phase], dft['distance']), 
                                                     vdims = ['Distance (degree)']
                                                     #label = phase
                                                    ) for phase in phaselist})
            # Make the phase curve thicker
            overlay_phase = overlay_phase * overlay_phase * overlay_phase * overlay_phase
        else:
            overlay_phase=None
    else:
        overlay_trace = [hv.Curve((dfm['time'], dfm[c]), vdims = ['Distance (degree)'], group='trace') for c in cols]
        if len(phaselist) > 0: 
            overlay_phase = [hv.Curve((dft[phase], dft['distance']), vdims = ['Distance (degree)'], group = 'phase') for phase in phaselist]
        else:
            overlay_phase=None
    
    return overlay_trace, overlay_phase, options

# read waveforms recorded by US Transportable Array
st = read('Files/waveforms/TA.*.LH*.mseed')

#overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [20,200], "Z", ["S", "pP", "P","PP"], "datashader")
overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [20,200], "Z", [], "datashader")
if overlay_phase != None:
    things_to_plot = overlay_trace * overlay_phase
else:
    things_to_plot = overlay_trace 
datashade(things_to_plot, cmap="gist_gray", cnorm='linear', aggregator=ds.count()).opts(options)

#overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [20,100], "T", ["S", "ScS"], "datashader", 1, 2)
overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [20,100], "T", [], "datashader", 1, 2)
if overlay_phase != None:
    things_to_plot = overlay_trace * overlay_phase
else:
    things_to_plot = overlay_trace 
datashade(things_to_plot, cmap="gist_gray", cnorm='linear', aggregator=ds.count()).opts(options)

TO DO#

Question 3.5 The record sections above show data from the first hour after the earthquake origin time on the Z and T components. Zoom into any distance range 10 degrees in width, and download the image as png from the options (top right). Draw a line with a red pen where you see the P wave (on the Z component) and S wave (on the T component). Find slope of this curve and report this as the speed of the P and S wave (in m/s). Are the slopes expected to vary with distance? If so, why? Hint: You can assume that 1 degree equals 111.1949 km on average. Attach the annotated image with the assignment.

Answer:

Part 2: Waveforms at different period bands#

Next, we will evaluate if we can glean more information by analyzing waveforms at other period bands. As you will see, there is a plethora of information in these seismograms at various period/frequency bands. Let us look at the vertical (Z) component from the US Transportable Array at two period bands (5-40 seconds vs 20-100 seconds).

overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [5,40], "Z", [], "datashader", 1, 2)
if overlay_phase != None:
    things_to_plot = overlay_trace * overlay_phase
else:
    things_to_plot = overlay_trace 
datashade(things_to_plot, cmap="gist_gray", cnorm='linear', aggregator=ds.count()).opts(options)

overlay_trace, overlay_phase, options = plottraces(cmt, st, 'Files/stations/', [20,100], "Z", [], "datashader", 1, 2)
if overlay_phase != None:
    things_to_plot = overlay_trace * overlay_phase
else:
    things_to_plot = overlay_trace 
datashade(things_to_plot, cmap="gist_gray", cnorm='linear', aggregator=ds.count()).opts(options)

TO DO#

Question 3.6 The record sections above show data from the first hour after the earthquake origin time on the Z component for data filtered at two period bands (5-40 seconds vs 20-100 seconds). Do you see any phases that are more prominent in one period band or the other? Annotate this phase by zooming into any distance range 10 degrees in width, and downloading the image as png from the options (top right). Attach the annotated image with the assignment.

Answer:

Part 3: Spectrum of Vibrations#

Mini-tutorial#

Primer in Fourier transforms

A Fourier transform (FT) is a mathematical transform that decomposes functions into frequency components, which are represented by the output of the transform as a function of frequency. Most commonly functions of time or space are transformed, which will output a function depending on temporal frequency or spatial frequency respectively.

%matplotlib inline

from scipy.fftpack import fft
import numpy as np
from importlib import reload
import os
from scipy import signal
import matplotlib.pyplot as plt

The Fourier transform takes a signal from the time into the frequency domain. One application of the Fourier transform is that we can recover the amplitudes and frequencies of a sampled signal. We will use the package numpy.fft. The underlying code for these functions is an f2c-translated and modified version of the FFTPACK routines. FFTPACK is a package of Fortran subprograms for the fast Fourier transform of periodic and other symmetric sequences. It includes complex, real, sine, cosine, and quarter-wave transforms.

For example, consider the signal 2 · cos(4 · 2πt) + 5 · sin(10 · 2πt) composed of a cosine with amplitude 2, frequency 4, and a sine with amplitude 5 and frequency 10. Using rfft of the package numpy.fft, the script below computes the discrete Fourier transform on the real array of samples via the efficient Fast Fourier Transform algorithm. We recover the amplitudes and corresponding frequencies of the components of our signal. With matplotlib we plot the spectrum.

f = lambda t: 2*np.cos(4*2*np.pi*t) + 5*np.sin(10*2*np.pi*t)
n = 256
r = np.arange(0.0,1.0,1.0/n)
s = f(r)
F = np.fft.rfft(s)
m = n/2
p = lambda z: (abs(np.real(z))/m,abs(np.imag(z))/m)
t = p(F)
tol = 1.0e-8
for i in range(0,len(t[0])):
    if(t[0][i] > tol):
        print(str(t[0][i]) + '*cos(' + str(i) + '*2*pi*t)')
    if(t[1][i] > tol):
        print(str(t[1][i]) + '*sin(' + str(i) + '*2*pi*t)')
fig = plt.figure()
plt.subplot(2, 1, 1)
plt.plot(abs(F)/m)
plt.subplot(2, 1, 2)
plt.plot(np.angle(F,deg=True))
plt.show()

2.0*cos(4*2*pi*t)
5.0*sin(10*2*pi*t)

To complete the illustration of time versus frequency plots (a normal signal versus it’s spectrum), here is a seismic signal (ground motion caused by an earthquake). The time series is a complicated looking signal. In fact, this signal can be described as a combination of a very large number of sinusoids, all with amplitudes as shown in the spectrum graph.

2004 Sumatra - Andaman Islands Earthquake USGS Link

The devastating M 9.1 earthquake off the west coast of northern Sumatra on December 26, 2004, occurred as the result of thrust faulting on the interface of the India plate and the Burma microplate. In a period of minutes, the faulting released elastic strains that had accumulated for centuries from ongoing subduction of the India plate beneath the overriding Burma microplate.

Since this earthquake was so large, it was able to ring the Earth like a bell, causing vibrational normal modes or standing waves. Let us look at the data from this earthquake recorded at Echery - Sainte Marie aux Mines, France.

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from obspy.clients.fdsn import Client
from obspy import UTCDateTime,read_inventory, read
from obspy.geodetics.base import gps2dist_azimuth
from gatspy.periodic import LombScargleFast

inventory_store = 'Files/ECH_inventory.xml'
path = Path(inventory_store)

if path.is_file():
    inv = read_inventory(inventory_store)
else:
    # Create two times which we will later use to download time series data.
    starttime = UTCDateTime(2004, 12, 26, 0, 58, 50)
    # Create a new time object by adding 10 days to the previous ones. Time
    # differences are in seconds.
    endtime = starttime + 86400 * 10
    client='IRIS'
    level='response'
    query ='G-ECH-*-LH*'

    c = Client(client)
    network,station,location,channel=query.split('-')
    inv = c.get_stations(network=network, station=station, location=location, channel=channel,
                     starttime=starttime, endtime=endtime,level=level)
    
    inv.write(inventory_store, format='STATIONXML')  

st = read('Files/Sumatra_2004_ECH.mseed')

nhours = 230
pre_filt=[1e-4,2e-4,3e-3,4e-3]

# read stat_lon, stat_lat
stat_lon = inv[0][0].longitude
stat_lat = inv[0][0].latitude

# read evt_lat, evt_lon
evt_lat = 3.09
evt_lon = 94.26

# compute back azimuth (check sign, later!):
dist, back_az, _ = gps2dist_azimuth(evt_lat, evt_lon,stat_lat, stat_lon)        

st2 = st.copy().slice(endtime=st[0].stats.starttime + 3600*nhours).remove_response(inventory=inv, output="ACC",pre_filt=pre_filt).rotate(method="NE->RT", inventory=inv,back_azimuth=back_az).taper(0.05)

# Remember the copy()!
tr = st2.select(component="Z")[0]


fmin = 0.0001
fmax = 0.001
N = 100000 # frequency resolution
df = (fmax - fmin) / N

time = np.linspace(start=0,stop=tr.stats.endtime-tr.stats.starttime,num=tr.stats.npts,endpoint=True)
model = LombScargleFast().fit(time, tr.data)
power = model.score_frequency_grid(fmin, df, N)
freqs = fmin + df * np.arange(N)

# convert to dataframe
data = dict()
data['Frequency [mHz]'] = freqs*1000 
data['Power'] = power
dfm = pd.DataFrame(data)


fig_title = 'Station = %s-%s, %s Channel' % (tr.stats.network,tr.stats.station,tr.stats.channel)
fig_1 = px.line(dfm, x='Frequency [mHz]', y='Power', 
                title=fig_title)
fig = go.Figure(data=fig_1,layout=fig_1.layout)
plotly.offline.iplot(fig)

TO DO#

Question 3.7 Report the approximate frequencies that you observe on the spectrum here. Identify atleast one normal mode and describe how it deforms the Earth.

Answer:

# LOOK AT TIDES
fmin = 0.001*1E-3
fmax = 0.03*1E-3

pre_filt=[1e-6,2e-6,3e-3,4e-3]
st2 = st.copy().slice(endtime=st[0].stats.starttime + 3600*nhours).remove_response(inventory=inv, 
                      output="ACC",pre_filt=pre_filt).rotate(method="NE->RT", 
                      inventory=inv,back_azimuth=back_az).taper(0.05)
tr = st2.select(component="Z")[0]
tr.plot()

N = 60 # frequency resolution
df = (fmax - fmin) / N

# Faster
time = np.linspace(start=0,stop=tr.stats.endtime-tr.stats.starttime,num=tr.stats.npts,endpoint=True)
model = LombScargleFast().fit(time, tr.data)
power = model.score_frequency_grid(fmin, df, N)
freqs = fmin + df * np.arange(N)

# convert to dataframe
data = dict()
data['Frequency [mHz]'] = freqs*1000 
data['Power'] = power
dfm = pd.DataFrame(data)


fig_title = 'Station = %s-%s, %s Channel' % (tr.stats.network,tr.stats.station,tr.stats.channel)
fig_1 = px.line(dfm, x='Frequency [mHz]', y='Power', 
                title=fig_title)
fig = go.Figure(data=fig_1,layout=fig_1.layout)
plotly.offline.iplot(fig)

TO DO#

Question 3.8 Seismometer are so sensitive that they also record vibrations due to tides in the solid Earth. An example is seen in the spectrum above. Report the two frequencies. Convert to minutes/hours and comment on the forces causing these longer period vibrations (compared to seismic waves).

Answer:

Fundamentals of Solid Earth Science

Seismology and Isostasy

Contents

Seismology and Isostasy#

Downloads from the Internet#

Part I: Mountain building and isostasy#

TO DO#

TO DO#

Part II: Seismic stations, seismographs, and seismograms#

Mini-tutorial#

TO DO#

Section A: Signal processing of vibrations at a single station#

Section B: Horizontal versus Vertical Motions#

Task I: At a single station#

TO DO#

Task II: Record Section of the USArray Transportable Array#

TO DO#

Part 2: Waveforms at different period bands#

TO DO#

Part 3: Spectrum of Vibrations#

Mini-tutorial#

TO DO#

TO DO#