Relationship between time difference and depth in 1D velocity model.#

Converting Ps time difference to depth is an useful method to estimate the depth of subsurface discontinuities corresponding to RF pulses. As the formula between Ps time difference \(T_{Ps}\) and depth (equivalent to radius \(R(i)\) in spherical coordinates) under a layered model, \(V_P\) and \(V_S\) influence the relationship.

\[T_{Ps}(n) = \sum_{i=1}^n\left(\sqrt{\left(\frac{R(i)}{V_S(i)}\right)^2-p_{Ps}^2}-\sqrt{\left(\frac{R(i)}{V_P(i)}\right)^2-p^2}\right)\frac{{\Delta}r}{R(i)}\]

In relevant studies, different velocity models are used to estimate the effects of velocity anormalies on depth of discontinuities.

Define a velocity model#

The velocity model should be prepared in a txt file with 3 columns (i.e., depth or thickness, Vp and Vs). Both depth or thickness are accepted for Seispy.

  • Using thickness of layers are popular in other toolbox, such as the CPS, Raysum and Hk.

  • The depth is defined in reference models and in the Taup toolkit.

from seispy.utils import DepModel
import matplotlib.pyplot as plt
import numpy as np

Customize a layered model#

Note

This is a new feature in version 1.2.19. Please update to dev branch for practice.

  • h is the thickness of layers. The last layer can be set to 0 as the half space.

  • vp and vs are velocities of corresponding layers.

h = np.array([20, 15, 175, 200, 250, 0])
vp = np.array([5800, 6500, 8175, 8665, 9864, 10923])/1000
vs = np.array([3460, 3850, 4500, 4783, 5398, 6089])/1000
dep_range = np.arange(800)
depmod_layer = DepModel.read_layer_model(dep_range, h, vp, vs)
depmod_layer.plot_model(show=False)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [2], line 5
      3 vs = np.array([3460, 3850, 4500, 4783, 5398, 6089])/1000
      4 dep_range = np.arange(800)
----> 5 depmod_layer = DepModel.read_layer_model(dep_range, h, vp, vs)
      6 depmod_layer.plot_model(show=False)

AttributeError: type object 'DepModel' has no attribute 'read_layer_model'

Use built-in reference model#

Seispy has imported some reference models. ak135 and iasp91 are available for argument velmod.

depmod_ak135 = DepModel(dep_range, velmod='ak135')
depmod_ak135.plot_model(show=False)
../_images/94e12d0edc9486438d66a1981c9772640e42337b6f19f8e5977394974252fac3.png

Calculate time difference#

The DepModel involve method for time difference calculation of conversion and multiple phases. Here, assuming the ray-parameter as 0.06 s/km, we calculate time differences with dep_range set before.

Note

The ray-parameters of P and Ps can be different. In following case, the ray-parameter of Ps is approximated as that of direct P wave.

from cProfile import label
from seispy.geo import skm2srad
rayp = skm2srad(0.06)
t1 = depmod_ak135.tpds(rayp, rayp) # Ps
tm1 = depmod_ak135.tpppds(rayp, rayp) # PpPs
tm2 = depmod_ak135.tpspds(rayp, rayp) # PsPs+PpSs

plt.plot(t1, depmod_ak135.depths, label='Ps')
plt.plot(tm1, depmod_ak135.depths, label='PpPs')
plt.plot(tm2, depmod_ak135.depths, label='PsPs+PpSs')
plt.legend()
plt.gca().invert_yaxis()
plt.xlabel('Time (s)')
plt.ylabel('Depth (km)')

# inset axes....
axins = plt.gca().inset_axes([0, 1.1, 0.47, 0.47])
axins.plot(t1, depmod_ak135.depths, label='Ps')
axins.plot(tm1, depmod_ak135.depths, label='PpPs')
axins.plot(tm2, depmod_ak135.depths, label='PsPs+PpSs')
# sub region of the original image
x1, x2, y1, y2 = 0, 30, 0, 150
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)
axins.invert_yaxis()
plt.gca().indicate_inset_zoom(axins, edgecolor="black")
(<matplotlib.patches.Rectangle at 0x17e1508b0>,
 (<matplotlib.patches.ConnectionPatch at 0x17e18ab20>,
  <matplotlib.patches.ConnectionPatch at 0x17e18aee0>,
  <matplotlib.patches.ConnectionPatch at 0x17e19a2e0>,
  <matplotlib.patches.ConnectionPatch at 0x17e19a580>))
../_images/7fe4c389d1573cc0b49296a2f6b505b53200da695375b4b52ba8625c6bee187b.png

Following case show relationship between time difference and depth with ray-parameter from 0.04 to 0.08 s/km.

rayp_list = np.arange(0.04, 0.09, 0.01)
tps = np.zeros([rayp_list.size, depmod_ak135.depths.size])
for i, raypps in enumerate(rayp_list):
    tps[i] = depmod_ak135.tpds(skm2srad(raypps), skm2srad(raypps)) 
    plt.plot(tps[i], depmod_ak135.depths, label='{:.2f}'.format(raypps))
plt.legend(title='Ray-parameter (s/km)')
plt.gca().invert_yaxis()
plt.xlabel('Time (s)')
plt.ylabel('Depth (km)')
/Users/xumijian/miniconda3_arm/lib/python3.9/site-packages/seispy/utils.py:156: RuntimeWarning: invalid value encountered in sqrt
  np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) *
Text(0, 0.5, 'Depth (km)')
../_images/9a0b8afe271c74f8e4cff3617970799d23008399e18712e6b0d2621c675b8c07.png