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)
A point source can be computed as a sum of all point sources weighted by their moment:
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)
Task: Compare the seismograms for three different representations of the 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.dt) * 50 finite_source.resample_sliprate(dt=finite_source.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()
# 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 = # 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 = # compute seismogram for finite source tr_finite = plt.plot(tr_simple.times(), tr_simple.data, label='simple') plt.plot(...) plt.plot(...) plt.legend() plt.xlim(0, tr_simple.times()[-1]) plt.show()