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

Note

This notebook can be downloaded as ps-depth.ipynb

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#

• 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)


### 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)


## 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 seispy.geo import skm2srad
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 0x108af6b00>,
(<matplotlib.patches.ConnectionPatch at 0x108b78820>,
<matplotlib.patches.ConnectionPatch at 0x108b78a60>,
<matplotlib.patches.ConnectionPatch at 0x108b78e80>,
<matplotlib.patches.ConnectionPatch at 0x108b79180>))


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):

/usr/local/miniconda/envs/docs/lib/python3.10/site-packages/seispy/utils.py:196: RuntimeWarning: invalid value encountered in sqrt

Text(0, 0.5, 'Depth (km)')