Ambiguity between hypocentral and station velocity when not using NLL grids

[krisvanneste]

I noticed that the _displacement_to_moment function in ssp_build_spectra calls ssp_util.get_vel to obtain vs_station.
However, when there is no NLL grid, ssp_util.get_vel just returns config.vs, which according to sourcespec.conf is the S-wave velocity close to the source. The _hypo_vel function in ssp_read_traces also calls ssp_util.get_vel and stores the same value in config.hypo.vs. So, when not using NLL grids, hypocentral and station velocities end up with the same value.

Would it be possible to define hypocentral and station velocities (P and S) separately in sourcespec.conf?

[claudiodsf]

I see the problem here, thank you.

So you propose to have two sets of config parameters in the INVERSION PARAMETERS section of config file (let me know if you have better defaults):

# P and S wave velocity close to the source (km/s)
vp_source = float(default=5.5)
vs_source = float(default=3.2)
# P and S wave velocity close to the stations (km/s)
vp_stations = float(default=5.5)
vs_stations = float(default=3.2)

Note that there is another set of velocities in the TIME WINDOW PARAMETERS section:

# P and S wave velocity (in km/s) for travel time calculation
# (if None, the global velocity model 'iasp91' is used)
vp_tt = float(default=None)
vs_tt = float(default=None)

But I guess that it's ok to leave those, since it's a sort of average velocity between source and stations.

[krisvanneste]

Claudio,

Yes, separate v[p/s]_source and v[p/s]_stations configuration parameters would be great. However, it means that older sourcespec.conf files will not work anymore... Another difficulty will be to update the ssp_util.get_vel function, which should return either v_source or v_stations. I suppose this can be done based on the depth value, because the hypocentral depth is known in config.hypo.depth.

Concerning vp_tt and vs_tt: Besides travel time calculation (in ssp_wave_arrival.add_arrivals_to_trace), these parameters are also used to compute the quality factor (in ssp_inversion._spec_inversion) and to set bounds on t_star (in ssp_data_types._Qo_to_t_star). In both cases, hypocentral velocities are used if they are not defined. Perhaps using the average of v_source and v_stations would be slightly better in that case.

Kris

[claudiodsf]

Yes, separate v[p/s]_source and v[p/s]_stations configuration parameters would be great. However, it means that older sourcespec.conf files will not work anymore...

It is not a goal of SourceSpec to keep backwards compatibility of sourcespec.conf: the config file can be updated via source_spec -U <config_file>, so this is a non-issue 😉

Another difficulty will be to update the ssp_util.get_vel function, which should return either v_source or v_stations. I suppose this can be done based on the depth value, because the hypocentral depth is known in config.hypo.depth.

Even better, I would compute, for the given arbitrary point, distance from station and distance from hypo and use the closest velocity.

Concerning vp_tt and vs_tt: Besides travel time calculation (in ssp_wave_arrival.add_arrivals_to_trace), these parameters are also used to compute the quality factor (in ssp_inversion._spec_inversion) and to set bounds on t_star (in ssp_data_types._Qo_to_t_star). In both cases, hypocentral velocities are used if they are not defined. Perhaps using the average of v_source and v_stations would be slightly better in that case.

Good catch!

I will make these modifications and report back here 😉

[krisvanneste]

Claudio,

OK, these modifications will probably impact the P_wave_type branch that I am preparing, but we will see.

Kris

[claudiodsf]

Don't worry, I will rebase before merging.

[claudiodsf]

I have just pushed 4 commits to the master branch which should fix this issue.

• concerning the automatic choice of vp[vs]_source with respect to vp[vs]_stations, I finally implemented your suggestion which was much easier 😉
• concerning quality factor computation, I now store travel time in trace and spectral objects and use it for the computation

Could you please check if the code works on your side?

P.S. For the moment, I assume the S-wave travel time. In your P_wave_type branch you should add a check...

[krisvanneste]

I just prepared the two other pull requests. I'm pretty sure the P_wave_type branch will conflict with your new modifications.
Unfortunately, I will not be able to test your changes and resolve the conflicts immediately, due to a meeting tomorrow.

[claudiodsf]

No problem: I also will be offline tomorrow (July 14 🇫🇷) and on Friday (bridge 😉).

[krisvanneste]

OK, enjoy your national holiday! Next week, it's our turn (July 21 + bridge as well)!

[krisvanneste]

Claudio,

I have closed the conflicting pull request and prepared a new one that is compatible with your new modifications. I did a test run with both P and SH waves and obtained exactly the same results as before. I also log the hypocentral and station velocities that are used in _displacement_to_moment, and they correspond to the ones I have configured. This confirms that your changes are working for me, but I must confess that I run sourcespec slightly differently than intended: I skip the parse_args, configure and read_traces stages, so if there would be problem there, I would not catch it.

[claudiodsf]

This confirms that your changes are working for me, but I must confess that I run sourcespec slightly differently than intended: I skip the parse_args, configure and read_traces stages, so if there would be problem there, I would not catch it.

Great! I will close this issue as resolved, then.

I'm very curious to know how you run SourceSpec (directly from Python?). Could you please share a use case?

[krisvanneste]

I'm very curious to know how you run SourceSpec (directly from Python?). Could you please share a use case?

Sure. Indeed, I can run sourcespec directly from Python.
I have written a generic interface between my own robspy library (extension of obspy) and sourcespec, with a function 'run_sourcespec' that takes an earthquake object, a list of seismograms, a station inventory, a phase pick dictionary and a number of optional arguments. You can have a look at it here: https://gitlab-as.oma.be/kris/robspy/-/blob/master/sourcespec.py

I also have a higher-level function that is part of a project in which we collect all available digital waveforms of earthquakes since 1985 in Belgium in order to construct a ground-motion database (BELSHAKE). In this function, I only have to pass the event ID and output folder (along with optional arguments) and only records with good quality for the selected wave type and component will be used. This code is not on gitlab yet, but I can add it if you are interested.

[claudiodsf]

Interesting!

Maybe one day (no rush!) we could include this Python interface into the code.

[krisvanneste]

Maybe one day (no rush!) we could include this Python interface into the code.

I think it would be possible to include some of it in SourceSpec directly. A first idea would be to add a function in source_spec.py that takes 'st' and 'config' arguments, and contains the code from the main section (minus parse_args, read_traces and ssp_exit) and perhaps returns proc_st, spec_st, specnoise_st, weight_st and sourcepar. This function could be called from the main section and perhaps it is also possible to add a higher-level function building further on that...

From my side, I will not be able to work on that before I return from holidays (mid-August).

[claudiodsf]

Sure, no need to rush on that.

(Working on "P_wave" PR right now ;-)