krischer / instaseis

Instant high-frequency seismograms from an AxiSEM database
http://instaseis.net
Other
49 stars 23 forks source link

P wave energy on the transverse component? #65

Closed rossmaguire closed 4 years ago

rossmaguire commented 5 years ago

Dear Instaseis team,

I have been using Instaseis along with the pre-computed Mars databases (http://instaseis.ethz.ch/marssynthetics/), and I noticed that there is a surprising amount of P wave energy on the transverse component. As an example, see the plot below for an Mw 5 event using Mars model DWTHot. The top panel shows both the radial and transverse plotted together, and the bottom shows a zoom-in on the transverse. The dashed vertical line is the TauP predicted P wave arrival. Data is filtered between 6 s and 50 s. In this case, the amplitude of the P wave arrival on the transverse is about 5 % of the radial, which seems too high. Any ideas why?

import instaseis
import numpy as np
import matplotlib.pyplot as plt

#read instaseis streams
db = instaseis.open_db('http://instaseis.ethz.ch/blindtest_5s/DWThot_5s/')
source = instaseis.Source(latitude=0.0,longitude=0.0,m_tt=-2.55e16,m_pp=2.55e16,m_tp=-3.04e16,depth_in_m=0)
receiver = instaseis.Receiver(latitude=0.0,longitude=20.0)
st = db.get_seismograms(source=source,receiver=receiver,kind='displacement',components='RTZ')
st.filter('bandpass',freqmin=1/50.0,freqmax=1/6.0)
tr_t = st.select(channel='LXT')[0]
tr_r = st.select(channel='LXR')[0]
time = np.linspace(0,tr_t.stats.npts*tr_t.stats.delta,tr_t.stats.npts)

#plot
xlim = [150,180]
ylim1 = [-3e-7, 3e-7]
ylim2 = [-2e-8, 2e-8]
P_arr = 159.58

fig,axes = plt.subplots(nrows=2,figsize=[8,5])
axes[0].plot(time,tr_t.data,label='transverse')
axes[0].plot(time,tr_r.data,label='radial')
axes[0].set_xlim(xlim)
axes[0].set_ylim(ylim1)
axes[0].axvline(P_arr,color='k',linestyle='dashed')
axes[0].legend(loc='upper left')
axes[0].set_ylabel('disp. (m)')
axes[0].ticklabel_format(axis='y',style='sci',scilimits=[-5,5])

axes[1].plot(time,tr_t.data)
axes[1].set_xlim(xlim)
axes[1].set_ylim(ylim2)
axes[1].axvline(P_arr,color='k',linestyle='dashed')
axes[1].set_xlabel('time (s)')
axes[1].set_ylabel('disp (m)')
axes[1].ticklabel_format(style='sci')

axes[0].set_title('http://instaseis.ethz.ch/blindtest_5s/DWThot_5s/, Mw 5 @ 20 degrees')
plt.show()

instaseis_dwthot_5s_rvst

Thanks!

krischer commented 5 years ago

I assigned this to @martinvandriel and @sstaehler who are involved in this database.

rossmaguire commented 5 years ago

Much appreciated!

martinvandriel commented 5 years ago

This is not a bug of this particular database but rather an unknown feature of theoretical seismology.

For comparison, replace the database with one of the syngine databases, e.g. db = instaseis.open_db('syngine://ak135f_2s')

image For a detailed reasoning, check Dahlen & Trump the section on Moment Tensor response: image As you can see, the tranverse component of a seismogram excited from a moment tensor source includes both toroidal and spheroidal motion, hence also P-waves. If you change the moment tensor to isotropic (e.g. only m_rr), these disappear because of the phi derivative in the displacement operator and you are back to what you assumed: spheroidal modes on radial and vertical, toroidal on transverse: image Cheers, Martin

rossmaguire commented 5 years ago

Hi Martin,

Thanks for the detailed response. It is very helpful.

If I understand your explanation of the implications of equation 10.60 from Dahlen and Tromp, the only way to completely separate the spheroidal and toroidal modes is to only include isotropic moment tensor terms (i.e. diagonal terms). Is this correct?

In my case, I am trying to generate synthetics with purely toroidial mode energy on the transverse component (i.e., I want zero energy from PSV related phases on my seismograms). Below I plotted synthetics from the DWThot_5s database for each moment tensor element (both isotropic and not). The source is at lon = 0, lat = 0, and the receiver is at lon = 20, lat = 0.

mars_mt1

Indeed if I use only m_rr, m_tt, or m_pp, the P wave is absent on the transverse, but all other energy is eliminated also. However, I think this is mostly due to the radiation pattern. For example, if I use a purely m_tt source, and move the receiver off of the equator (in this case at lon = 14, lat = 14), the P arrival shows up on the transverse. Again, it is about 5 % of the amplitude of the radial (see the figure below).

mars_mtt

At first glance, it looks like using an m_rt or m_tp source has no energy arriving at P but if I zoom in closer, it is still there. The plot below shows the synthetics for each moment tensor element, but this time zoomed in on the expected P wave arrival.

mars_mt2

I have also tried to use a single force in the theta direction, and the P wave still appears on the transverse.

mars_theta_force

In short, there doesn't seem to be any source I can use to make synthetics that include the SH wavefield, but don't include any P wave energy. Do you have any advice?

PS, Sorry for the deluge of plots.

martinvandriel commented 5 years ago

Hi,

the only way to completely separate the spheroidal and toroidal modes

to be precise: they are of course decoupled (if this is what you mean with separate), but both generate transverse and radial motion for non-isotropic sources. In a normal mode framework it would of course be easy to just sum the toroidal modes, but in instaseis that is not possible.

The source you are looking (isotropic, SH waves only) can be thought of as a rotational source. I do recall having seen images of such a device for exploration seismics, but can't find it just now. It was also used for numerical solutions by [1, 2] and in fact requires an antisymmetric moment tensor, but the symmetry of the moment tensor is a consequence of the preservation of angular momentum in earthquakes, so this is not observed in nature.

Instaseis hence does not have such a source built in, but you could mimic it by composing the moment tensor of the individual forces offset by a finite difference (say a few meters). The pure SH source than looks like a swastika from top.

cheers, Martin

[1] Jahnke, Gunnar, Michael S. Thorne, Alain Cochard, and Heiner Igel. 2008. “Global SH-Wave Propagation Using a Parallel Axisymmetric Spherical Finite-Difference Scheme: Application to Whole Mantle Scattering.” Geophysical Journal International 173 (3): 815–26. doi:10.1111/j.1365-246X.2008.03744.x.

[2] Igel, Heiner, and Michael Weber. 1995. “SH-Wave Propagation in the Whole Mantle Using High-Order Finite Differences.” Geophysical Research Letters 22 (6): 731–34.

rossmaguire commented 5 years ago

I understand. Thanks so much for taking the time to help!