Advanced Exercise 2: Finite Source Effects¶

For earthquakes with Magnitude up to about 5 recorded at teleseismic distances, approximating the fault by a point source is a reasonable approach. However, for larger earthquakes with longer rupture duration this approximation is not valid anymore. In this exercise, you will compare the point source approximation with finite source solutions to understand its limitations.

For three of the earthquakes we use in this tutorial, USGS provides finite fault solutions: the recent event in Nepal, the largest aftershock and the one in Chile. This is the fault solution and slip as a function of time for the Nepal M7.9 event:

Fault representation (image: USGS)

Source Time Function (image: USGS)

Basic lines to set up the notebook and some paths.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
import obspy
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10, 8)

Import Instaseis and open the database:

import instaseis
db = instaseis.open_db("data/database")

In Instaseis, a finite fault is represented as set of point sources, where each point source represents one of the fault patches with individual source time function. This functionality is provided by the instaseis.FiniteSource object (see Documentation). It can be initialized in two ways: from a list of point sources, or more conveniently by reading *.param files provided by USGS or standard rupture format (*.srf) files (these can also be used in the GUI).

finite_source = instaseis.FiniteSource.from_usgs_param_file(
    'data/events/finite_source/FINITE_SOURCE_2015_05_12__Mw_7_2_Nepal.param')
print(finite_source)

Instaseis Finite Source:
	Moment Magnitude     : 7.28
	Scalar Moment        :   1.04e+20 Nm
	#Point Sources       : 525
	Rupture Duration     :   35.8 s
	Time Shift           :    0.0 s
	Min Depth            : 8695.7 m
	Max Depth            : 20897.6 m
	Hypocenter Depth     : 20897.6 m
	Min Latitude         :   27.1 deg
	Max Latitude         :   28.3 deg
	Hypocenter Latitude  :   27.7 deg
	Min Longitude        :   85.4 deg
	Max Longitude        :   86.9 deg
	Hypocenter Longitude :   86.2 deg

A point source can be computed as a sum of all point sources weighted by their moment:

finite_source.compute_centroid()
print(finite_source.CMT)

Instaseis Source:
	Origin Time      : 1970-01-01T00:00:00.000000Z
	Longitude        :   86.2 deg
	Latitude         :   27.7 deg
	Depth            : 1.5e+01 km km
	Moment Magnitude :   7.26
	Scalar Moment    :   9.76e+19 Nm
	Mrr              :   2.63e+19 Nm
	Mtt              :  -2.46e+19 Nm
	Mpp              :  -1.70e+18 Nm
	Mrt              :   9.33e+19 Nm
	Mrp              :  -8.13e+18 Nm
	Mtp              :   9.86e+18 Nm

The hypo- and epicenter can be found as the fault patch that ruptures first:

finite_source.find_hypocenter()
print('hypocenter latitude:', finite_source.hypocenter_latitude,
      'longitude:', finite_source.hypocenter_longitude,
      'depth:', finite_source.hypocenter_depth_in_m / 1e3)

hypocenter latitude: 27.664442081754583 longitude: 86.1627 depth: 14.7966

Task: Compare the seismograms for three different representations of the source:

A point source with simple gaussian source time function (using CMTSOLUTION or quakeml files),
the CMT solution using the more complex source time function provided by finite_source.CMT
the full finite source solution using the FiniteSource object and db.get_seismograms_finite_source()

Note: First, you have to adapt the sampling of the source time functions in the finite source to the database, which works like this:

# reloading finite source here to be sure to have a clean source time function
finite_source = instaseis.FiniteSource.from_usgs_param_file(
    'data/events/finite_source/FINITE_SOURCE_2015_05_12__Mw_7_2_Nepal.param')

# prepare the source time functions to be at the same sampling as the database
# first use enough samples such that the lowpassed stf will still be correctly represented
nsamp = int(db.info.period / finite_source[0].dt) * 50
finite_source.resample_sliprate(dt=finite_source[0].dt, nsamp=nsamp)
# lowpass to avoid aliasing
finite_source.lp_sliprate(freq=1.0/db.info.period)
# finally resample to the sampling as the database
finite_source.resample_sliprate(dt=db.info.dt, nsamp=db.info.npts)

finite_source.compute_centroid()

Solution

# load receivers from stations xml file
receivers = instaseis.Receiver.parse('data/stations/all_stations.xml')

simple_source = instaseis.Source.parse(
    'data/events/quakeml/GCMT_2015_04_25__Mw_7_9.xml')

# compute seismogram with CMT solution and no simple source time function (gaussian):
tr_simple = db.get_seismograms(
    simple_source, receivers[0], components=('Z'), dt=1.0)[0]

# compute seismogram with CMT solution and source time function computed as the 
# sum of all source time functions in the finite source (reconvolve_stf=True):
tr_cmt = db.get_seismograms(
    finite_source.CMT, receivers[0], components=('Z'),
    dt=1.0, reconvolve_stf=True, remove_source_shift=False)[0]

# compute seismogram for finite source
tr_finite = db.get_seismograms_finite_source(
    finite_source, receivers[0], components=('Z'), dt=1.0)[0]

plt.plot(tr_simple.times(), tr_simple.data, label='simple')
plt.plot(tr_cmt.times(), tr_cmt.data, label='cmt')
plt.plot(tr_finite.times(), tr_finite.data, label='finite')
 
plt.legend()
plt.xlim(0, tr_simple.times()[-1])

plt.show()