Travel time and takeoff angle from 1D velocity model

[claudiodsf]

This PR replaces PR #78 by implementing the changes in SourceSpec V1 instead of V2. See also issue #57.

I’ll port everything to V2 once this PR is merged, which will involve a fairly large (and much-needed) rebase.

@krisvanneste, I edited ssp_velocity_model.py to better follow SourceSpec coding conventions and to handle some errors more cleanly by raising exceptions. You can see the details by diffing this version against the previous one.

Note that I haven’t been able to test the code yet, since the minimize_scalar() function is currently missing: could you please commit this function to this branch?

[krisvanneste]

Claudio, thanks.

I just had a look at the missing minimize_scalar function, but this should be present in scipy.optimize, in at least versions 1.6.2 through 1.13.0. Which version do you use?

[claudiodsf]

Ah ok, thanks. There was no import in your file. I'm adding this and pushing a new version 😉

[claudiodsf]

Force pushed! Do you mind checking if the code works as intended?

[krisvanneste]

Ah ok, thanks. There was no import in your file. I'm adding this and pushing a new version 😉

I didn't notice, but I must have deleted this accidentally...

[krisvanneste]

Force pushed! Do you mind checking if the code works as intended?

I will try tomorrow or Wednesday.

[claudiodsf]

Hi Kris, so, I could start testing, and I found a basic issue.

You get the velocity model from vp_source and vs_source, but this is not consistent with what is explained in the config file:

# TIME WINDOW PARAMETERS --------
# P and S wave velocity (in km/s) for travel time calculation
# (if None, the global velocity model 'iasp91' is used)
# Theoretical P or S arrival times are used when a manual P or S pick is not
# available, or when the manual P or S pick is too different from the
# theoretical arrival (see 'p_arrival_tolerance' and 's_arrival_tolerance'
# below).
vp_tt = None
vs_tt = None

So, I think we should extend vp_tt and vs_tt to accept a list of velocities. We should also add something like layer_top_depths_tt to indicate the layer depths for travel time calculation.

Finally, we should modify the behaviour for when vp_source and vs_source are set to None. Currently:

# Layer top depths (km, positive down), for layered models (see below)
#  Note: generally, the first layer top depth should be 0 or a negative value
layer_top_depths = force_list(default=None)
# P and S wave velocity close to the source (km/s)
# It can be a single value or a list of values (layered model)
# Set to None to use velocity from the global Earth model 'iasp91'
#   Note: specifying a layered model is useful when the same config file is
#   used for several SourceSpec runs with sources at different depths
vp_source = force_list(default=5.5)
vs_source = force_list(default=3.2)

I propose to replace

# Set to None to use velocity from the global Earth model 'iasp91'

with

# Set to None to use 'vp_tt' and 'vs_tt' if available (see above); if not, fall back to velocities from the global Earth model 'iasp91'.

What do you think?

[claudiodsf]

Ok, I pushed a new commit implementing my proposal above.

The relevant part in config file is:

# TIME WINDOW PARAMETERS --------
# Layer top depths (km, positive down), for travel time calculation.
# Only needed if more than one velocity value is provided for 'vp_tt' and
# 'vs_tt' below (layered model).
#   Note: generally, the first layer top depth should be 0 or a negative value
layer_top_depths_tt = force_list(default=None)
# P and S wave velocity (in km/s) for travel time calculation.
# These can be specified either as a single value or as a list of values
# (layered model; see 'layer_top_depths_tt' above).
# If set to None, the global velocity model "iasp91" is used.
#   Note: Theoretical P or S arrival times are used when a manual P or S pick is
#   not available, or when the manual P or S pick is too different from the
#   theoretical arrival (see 'p_arrival_tolerance' and 's_arrival_tolerance'
#   below).
vp_tt = force_list(default=None)
vs_tt = force_list(default=None)

Let me know if that's ok for you.

Also note that, in this new implementation, the function _wave_arrival_vel() (Travel time and takeoff angle using a constant velocity) is never used, since it is superseded by _wave_arrival_layer_model(). Ok to remove it?

[krisvanneste]

Claudio, sorry, but I didn't have the time to look at this today. I will try tomorrow.

[claudiodsf]

Claudio, sorry, but I didn't have the time to look at this today. I will try tomorrow.

Sure! No problem.

I added an extra commit implementing the fallback to vp_tt and vs_tt in case vp_source and vs_source are not defined.

[krisvanneste]

Claudio, I ran a quick test and it seems to work: if I don't pass any phase picks, theoretical travel times are calculated from the 1D velocity model and these appear to be reasonable.

Concerning the configuration file, the current mechanism is working, but maybe it would be more logical if there is only 1 way to define a layered velocity model? Now both v[p/s]_source and v[p/s]_tt can be lists, which may be confusing. I think it would be more clean if there would be a general velocity model, and v[p/s]_source, v[p/s]_tt and v[p/s]_station are derived from that if they are not specified individually. Maybe grouping all velocity-related configuration options would make this more obvious. This is just an idea, not a big deal for me.

[claudiodsf]

Thanks, Kris, for the testing and suggestions.

I just pushed a commit which consolidates the earth model parameters in the config file.

The new section is as follows:

# EARTH MODEL PARAMETERS --------
# General Earth model parameters.
# The parameters are defined in terms of vp (km/s), vs (km/s)
# and density (rho, kg/m3).
# These values can be provided either as single values or as layered models
# (list of values). In the latter case, the top depths of each layer
# (in km, positive down) must be provided as well.
#
# Earth model parameters are used to compute theoretical P and S
# arrival times as well as coefficients for converting station displacements
# to seismic moment.
#
# If any of the parameters below is set to None, values from the global
# Earth model "iasp91" are used.
#
# Example of layered model:
#   vp = 5.30, 5.60, 6.20, 6.90, 7.40, 7.70, 7.90, 8.10, 8.30
#   vs = 3.01, 3.18, 3.52, 3.92, 4.20, 4.37, 4.49, 4.60, 4.72
#   rho = 2520, 2610, 2780, 2970, 3120, 3200, 3260, 3320, 3370
#   layer_top_depths = 0.0, 4.0, 9.0, 14.0, 19.0, 24.0, 33.0, 49.0, 66.0
vp = force_list(default=5.5)
vs = force_list(default=3.2)
rho = force_list(default=2500)
layer_top_depths = force_list(default=None)
# Source-specific parameters.
# You can specify different Earth model parameters near the source.
# These values can be provided either as single values or as layered models
# (list of values). In the latter case, the top depths of each layer
# (in km, positive down) must be provided as well.
# If any of the parameters below is set to None, then general Earth model
# values (defined above) are used.
vp_source = force_list(default=None)
vs_source = force_list(default=None)
rho_source = force_list(default=None)
layer_top_depths_source = force_list(default=None)
# Station-specific parameters.
# You can specify different Earth model parameters near the stations.
# If any of the parameters below is set to None, then general Earth model
# values (defined above) are used.
# Note:
#   Station-specific parameters cannot be layered models.
vp_stations = float(min=0, default=None)
vs_stations = float(min=0, default=None)
rho_stations = float(min=0, default=None)
# NonLinLoc models.
# As an alternative, NonLinLoc valocity and travel time models can be
# specified.
# NLL_time_dir: directory containing NonLinLoc travel time grids.
# Theoretical travel times will be read from these grids, instead of being
# computed from the vp and vs values defined above.
# Note:
#   Reading NonLinLoc grids takes time. For simple 1D models, you
#   can speed up considerably the process using a generic station
#   named "DEFAULT". The coordinates of this default station are not important,
#   since they will be superseded by each station's coordinates.
NLL_time_dir = string(default=None)
# NLL_model_dir: directory containing a NonLinLoc velocity model.
# This model will be used to read vp, vs values near the source and stations,
# instead of using the values defined above.
# Important note:
#   NonLinLoc velocity models do not contain density values. Therefore,
#   when using a NonLinLoc velocity model, density values must be provided
#   through the parameters defined above (rho, rho_source, rho_stations).
NLL_model_dir = string(default=None)
# -------- EARTH MODEL PARAMETERS

Note that all the _tt parameters have disappeared and replaced by simple vp, vs, rho, layer_top_depths.

I ran some quick test, but this needs to be tested more thoroughly.

Let me know what you think!

[krisvanneste]

This looks great to me and quite logical. If I understand correctly, it is still possible to specify just 1 vp/vs value to calculate travel times. In that case, do you use _wave_arrival_vel() or does it also work with wave_arrival_layer_model()? I didn't test the latter case.

[claudiodsf]

If I understand correctly, it is still possible to specify just 1 vp/vs value to calculate travel times. In that case, do you use _wave_arrival_vel() or does it also work with wave_arrival_layer_model()? I didn't test the latter case.

Your code technically works for scalar values, but I didn't check the results. If they're compatible with _wave_arrival_vel(), we can safely remove this latter 😉

[krisvanneste]

If they're compatible with _wave_arrival_vel(), we can safely remove this latter 😉

OK, but we need to test this.

[claudiodsf]

If they're compatible with _wave_arrival_vel(), we can safely remove this latter 😉

OK, but we need to test this.

It looks like the results are the same 😉 :

P dist: 10km, z: 0km tt_homog=1.82s, tt_layer=1.82s, tt_diff=0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
P dist: 10km, z: 10km tt_homog=2.57s, tt_layer=2.57s, tt_diff=0.00s, a_homog=135.00°, a_layer=135.00°, a_diff=0.00°
P dist: 10km, z: 30km tt_homog=5.75s, tt_layer=5.75s, tt_diff=0.00s, a_homog=161.57°, a_layer=161.57°, a_diff=0.00°
P dist: 50km, z: 0km tt_homog=9.09s, tt_layer=9.09s, tt_diff=0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
P dist: 50km, z: 10km tt_homog=9.27s, tt_layer=9.27s, tt_diff=0.00s, a_homog=101.31°, a_layer=101.31°, a_diff=-0.00°
P dist: 50km, z: 30km tt_homog=10.60s, tt_layer=10.60s, tt_diff=0.00s, a_homog=120.96°, a_layer=120.96°, a_diff=0.00°
P dist: 100km, z: 0km tt_homog=18.18s, tt_layer=18.18s, tt_diff=0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
P dist: 100km, z: 10km tt_homog=18.27s, tt_layer=18.27s, tt_diff=0.00s, a_homog=95.71°, a_layer=95.71°, a_diff=0.00°
P dist: 100km, z: 30km tt_homog=18.98s, tt_layer=18.98s, tt_diff=0.00s, a_homog=106.70°, a_layer=106.70°, a_diff=0.00°
S dist: 10km, z: 0km tt_homog=3.12s, tt_layer=3.12s, tt_diff=-0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
S dist: 10km, z: 10km tt_homog=4.42s, tt_layer=4.42s, tt_diff=-0.00s, a_homog=135.00°, a_layer=135.00°, a_diff=0.00°
S dist: 10km, z: 30km tt_homog=9.88s, tt_layer=9.88s, tt_diff=-0.00s, a_homog=161.57°, a_layer=161.57°, a_diff=0.00°
S dist: 50km, z: 0km tt_homog=15.62s, tt_layer=15.62s, tt_diff=-0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
S dist: 50km, z: 10km tt_homog=15.93s, tt_layer=15.93s, tt_diff=-0.00s, a_homog=101.31°, a_layer=101.31°, a_diff=-0.00°
S dist: 50km, z: 30km tt_homog=18.22s, tt_layer=18.22s, tt_diff=-0.00s, a_homog=120.96°, a_layer=120.96°, a_diff=0.00°
S dist: 100km, z: 0km tt_homog=31.25s, tt_layer=31.25s, tt_diff=-0.00s, a_homog=90.00°, a_layer=90.00°, a_diff=0.00°
S dist: 100km, z: 10km tt_homog=31.41s, tt_layer=31.41s, tt_diff=-0.00s, a_homog=95.71°, a_layer=95.71°, a_diff=0.00°
S dist: 100km, z: 30km tt_homog=32.63s, tt_layer=32.63s, tt_diff=-0.00s, a_homog=106.70°, a_layer=106.70°, a_diff=0.00°

The script is here test_wave_arrival.py, if you want to play around a little more.

Note that for this constant velocity model, I provided a layer top depth of -999 km to your function, since a top depth of 0 km generates an error for an event at 0 km of depth (need to fix that in ssp_wave_arrival()).

[claudiodsf]

Note that for this constant velocity model, I provided a layer top depth of -999 km to your function, since a top depth of 0 km generates an error for an event at 0 km of depth (need to fix that in ssp_wave_arrival()).

This makes me wonder if we don't systematically assign a very large negative top depth to the first layer in a 1d model, to avoid errors for zero-depth or negative-depth earthquakes (everything can happen !)

[krisvanneste]

Note that for this constant velocity model, I provided a layer top depth of -999 km to your function, since a top depth of 0 km generates an error for an event at 0 km of depth (need to fix that in ssp_wave_arrival()).

This makes me wonder if we don't systematically assign a very large negative top depth to the first layer in a 1d model, to avoid errors for zero-depth or negative-depth earthquakes (everything can happen !)

Claudio,

Your test is pretty convincing.

However, the fix is only necessary for the constant velocity model. If more than 1 layer is provided, the layered model also works for zero or negative depths. If you would specify the top depth as -999 km, the travel times will be very long!

[krisvanneste]

One reason for keeping the original _wave_arrival_vel function is that it is significantly faster:

timeit _wave_arrival_vel(trace, vel[phase][0])
1.53 µs ± 64.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

timeit _wave_arrival_layer_model(trace, phase, vel)
146 µs ± 1.28 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

[claudiodsf]

Ok, so I kept both. Just force-pushed a new commit. Needs testing ;-)

[krisvanneste]

Needs testing ;-)

I obtain exactly the same results as in my previous test.

[claudiodsf]

I obtain exactly the same results as in my previous test.

Great! Let me finalize a couple of things (config file migration and changelog) and then I will merge this.

[claudiodsf]

This is complete from my point of view.

Please let me know if you see other things to change/improve!