In this example we generate P receiver functions for a model that includes a lower-crustal anisotropic layer. This example follows that of Figure 2 in Porter et al. (2011), which uses the Raysum software developed by Frederiksen and Bostock (2000)
Start by importing the necessary modules
import numpy as np
from obspy.core import Stream
from telewavesim import utils as ut
from telewavesim import wiggle as wg
Select the model file:
modfile = '../models/model_Porter2011.txt'
Select the type of incident wave - options are 'P'
, 'SV'
, 'SH'
, or 'Si'
, which is an isotropic S-wave source
Danger! Using 'SH' will not work properly for modeling receiver functions as the code will think you want plane-wave displacements (see below). Do not use 'SH' if you want S-wave receiver functions.
wvtype = 'P'
Next we use variables to define the desired time series.
npts = 3000 # Number of samples
dt = 0.01 # Sample distance in seconds
Now specify the parameters of the incident wavefield in terms of a horizontal slowness and back-azimuth. In this example the slowness won't change, so we can pass it as a global variable now. The back-azimuth will range from 0 to 360 degrees with a 10-degree increment, so we define a np.ndarray
for this variable and do not yet pass it as a global variable.
slow = 0.06 # Horizontal slowness (or ray parameter) in s/km
baz = np.arange(0., 360., 10.)
Read the model parameters and return a Model object. Up to here, the steps could have been performed in no particular order, except the name of the file that needs to be defined before the call to read_model()
model = ut.read_model(modfile)
As we need to loop through back-azimuth values, we will initialize empty Stream
objects to store the traces from the output of the main routine.
trR = Stream(); trT = Stream()
Now the main loop over back-azimuths where all calculations are done - self explanatory
Remember that the obs
boolean variable defaults to False
, so if you want to change to True
, either explicitely set them prior to this step or use the following call to ut.run_plane()
with argument obs=True
. Here we are not simulating OBS seismograms, so we don't need to specify anything.
# Loop over range of data
for bb in baz:
# Calculate the plane waves seismograms
trxyz = ut.run_plane(model, slow, npts, dt, bb, wvtype=wvtype, obs=False)
# Then the transfer functions in Z-R-T coordinate system
tfs = ut.tf_from_xyz(trxyz, pvh=False)
# Append to streams
trR.append(tfs[0]); trT.append(tfs[1])
The result of the previous loop is a set of transfer functions. To get receiver functions, we simply filter the Stream
objects using some frequency corners
# Set frequency corners in Hz
f1 = 0.01
f2 = 1.0
# Filter to get wave-like traces
trR.filter('bandpass',freqmin=f1, freqmax=f2, corners=2, zerophase=True)
trT.filter('bandpass',freqmin=f1, freqmax=f2, corners=2, zerophase=True)
36 Trace(s) in Stream: ... | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:29.990000Z | 100.0 Hz, 3000 samples ... (34 other traces) ... ... | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:29.990000Z | 100.0 Hz, 3000 samples [Use "print(Stream.__str__(extended=True))" to print all Traces]
Now plot the result as wiggles, using the format displayed in the paper by Porter et al. (2011). We also need to define 'stacked traces' that represent the average of all recceiver functions to be displayed.
# Stack over all traces
trR_stack, trT_stack = ut.stack_all(trR, trT, pws=True)
# Plot as wiggles
wg.rf_wiggles_baz(trR, trT, trR_stack, trT_stack, 'test', btyp='baz',
scale=1.e3, tmin=-5., tmax=8., save=False, ftitle='porter2011',
wvtype='P')
Stacking ALL traces in streams Plotting Wiggles by baz
And Voilà!
Now try the same example but setting wvtype = 'SV'
to get S receiver functions for the same model. Such an example would correspond to a core-refracted shear wave (such as SKS) with no incident transverse component. In this case the receiver functions are unstable (spectral division by zero-valued SH component) and are not computed for the transverse component. Setting wvtype = 'Si'
will produce a transverse component receiver function, which is more realistic for incident S waves propagating through the mantle only. Be careful with the slowness values!!!
Did you notice the boolean pvh=False
in the code above? It sets whether or not the seismograms are rotated to the P-SV-SH wave modes. Setting it to True
will essentially make the zero-lag signal on the radial. component disappear, since we know the exact value of the seismic velocities at the surface. This may not always be true so use with caution when comparing with real data!