Obspy: Tau-p underside reflections in lower mantle.

[JackWalpole]

Dear obspy team,

Would it be possible to improve the taup module so that it can get the ray paths for underside reflections in the lower mantle? Currently I am not able to get underside reflections for any interface deeper than the 660. When deeper reflectors are requested (e.g. "S^1000S") the 660 path is returned. Being able to get lower mantle underside reflector paths would be a useful feature.

Cheers,

Jack

-System:

• ObsPy 1.0.2, Python 2.7, Platform OSX

• Installed with anaconda

  • Working Example:
    (Note that the "S^1000S" path appears to reflect from the 660 km discontinuity -- not the 1000 km requested using the nomenclature.)

from obspy.taup import TauPyModel
model = TauPyModel(model="iasp91")
arrivals = model.get_ray_paths(
... source_depth_in_km=50,
... distance_in_degree=90,
... phase_list=["SKS","SS","S^1000S"])
arrivals.plot()

[megies]

CC @crotwell

[ThomasLecocq]

I think it works if you begin by discretising the base earth model in 1km values in a new taup_model, then "build it" and run your custom model.

[seisman]

According to "Phase Naming in TauP":

The numbers given need not be the actual depth, the closest depth corresponding to a discontinuity in the model will be used. For example, if the time for P410s is requested in a model where the discontinuity was really located at 406.7 kilometers depth, the time returned would actually be for P406.7s.

In your case, since there is no such 1000 discontinuity, the closest 660 km discontinuity is used.

To get exactly S^1000S phase, you can make a copy of iasp91 reference model and manually add a discontinuity at the depth of 1000 km. For example, the original iasp91 model is:

  908.500   11.3140    6.3199    4.5173
  958.000   11.3960    6.3546    4.5459 
  1007.500   11.4761    6.3883    4.5744                           
  1057.000   11.5543    6.4211    4.6028

You can add a discontinuity at 1000 km (note the slightly different Vp and Vs at 1000 km):

   908.500   11.3140    6.3199    4.5173        
   958.000   11.3960    6.3546    4.5459   
   1000.000   11.4531    6.3763    4.5744 
   1000.000   11.4561    6.3783    4.5744 
   1007.500   11.4761    6.3883    4.5744 
   1057.000   11.5543    6.4211    4.6028

Then, you can use this modified model to get the S^1000S phase:

from obspy.taup.taup_create import build_taup_model
from obspy.taup import TauPyModel
build_taup_model("myiasp91.tvel", output_folder=".")
model = TauPyModel(model="myiasp91.npz")
arrivals = model.get_ray_paths(source_depth_in_km=50,  
                               distance_in_degree=90, 
                               phase_list=["SKS","SS","S^1000S"])
arrivals.plot()

[ThomasLecocq]

I think this:

model = TauPyModel(model="myiasp91.npz")

should be

model = TauPyModel(model="myiasp91")

without the npz

[seisman]

model = TauPyModel(model="myiasp91.npz")

works.

And

model = TauPyModel(model="myiasp91")

fails.

I'm using obspy 1.0.3 and Python 3.6.0 (anaconda 4.3.1).

[seisman]

I think it's better if both can be supported, and even better if TauPyModel can build the npz file automatically if the model is not builtin.

[seisman]

@ThomasLecocq I use output_folder="." when calling build_taup_model, so the output npz is in current directory. This explains why model = TauPyModel(model="myiasp91.npz") works and the other fails. I think it's related to #1714.

[crotwell]

Just as far as the phase naming is concerned, you have it figured out. Let me know if more questions.

I have thought about adding something in TauP so that the discontinuity had to be within X km of the phase name, but never got around to implementing the check or deciding of the right X. Perhaps worth doing, assuming a reasonable value of X could be found.

[JackWalpole]

Thanks for your answers and especially for the example provided by @seisman

[JackWalpole]

Here's a little bit of rough and ready code that is useful for automatically adding a "discontinuity" in the velocity model text file.

import numpy as np
def add_a_discontinuity_to_tvel_file(infile,outfile,depth):
    ### a function to add a discontinuity to a tvel model at arbitrary depth
    # read in the model from text file (assuming 2 header lines)
    a = np.genfromtxt(infile,dtype=None,skip_header=2)
    # interpolate the values at required depth
    vpz = np.interp([depth],a[:,0],a[:,1])
    vsz = np.interp([depth],a[:,0],a[:,2])
    vsrho = np.interp([depth],a[:,0],a[:,3])
    # check that layer is not already present at this depth
    if (not np.in1d([depth],a[:,0]).all()):
        # if not then add it
        rowdeeper = np.argmax(a[:,0]>depth)
        a = np.insert(a,rowdeeper,[[depth,vpz,vsz,vsrho]],axis=0)
    # make the discontinuity
    eta=0.00001 # small number to be added to values    
    rowdeeper = np.argmax(a[:,0]>depth)
    a = np.insert(a,rowdeeper,[[depth,vpz+eta,vsz+eta,vsrho+eta]],axis=0)
    # save to outfile
    np.savetxt(outfile,a,fmt='%12.6f',header='dummyline\ndummyline')

[krischer]

I think it's better if both can be supported, and even better if TauPyModel can build the npz file automatically if the model is not builtin.

This is a conscious omission in ObsPy - I don't think we really want people accidentally do that in a loop or something as that is by far the slowest part of everything. Model building should IMHO be an explicit function so people are aware that this is not something they should constantly do.

I'm closing this issue at it is IMHO solved by @JackWalpole's script. This issue will exist so people should be able to find it.