1 Introduction

One of the main functions involved with global monitoring of nuclear explosions is the development of Earth models for predicting seismic travel times. The ability to accurately predict travel times allows for more accurate calculation of seismic event locations. The development of higher-dimensional (i.e., 2D, 3D) global and regional Earth models for the crust and upper mantle has been an on-going process, with the goal for some to better image Earth structure using standard body-wave, surface wave, or joint-inversion techniques (e.g., Amaru 2007; Burdick et al. 2008; Chang et al. 2010; Li et al. 2008; Simmons et al. 2010; Steck et al. 2011; Syracuse et al. 2017, among others), while for others the primary focus is to better predict seismic travel times for event location (e.g., Ballard et al. 2016a; Myers et al. 2010; Phillips et al. 2007; Simmons et al. 2012).

With the ease of online access to many seismic networks and bulletins including the Global Seismic Network (GSN) and the International Seismological Centre (ISC), as well as the utilization of existing global, regional and local networks (e.g., collections by the United States Geological Survey and the International Monitoring System (IMS) [public access possible via virtual Data Exploitation Centre (vDEC) [https://www.ctbto.org/specials/vdec (Vaidya et al. 2009)]), the need and interest to quickly and accurately locate all seismic events has pushed magnitude thresholds lower and lower. As these thresholds decrease, the reliance on more regional-distance seismic stations (i.e., ≲ 15–18°) is increasing in order to have a sufficient number of associated seismic phases with which to locate the smaller magnitude events.

Most seismic monitoring agencies rely on fast travel time calculations based on distance-dependent, one-dimensional (1D) Earth models (varying with depth) in order to locate seismic events in near real-time. However, the complexities of the crust and upper mantle are problematic for seismic event locations when including regional seismic phases (i.e., Pn, Pg, Sn, Lg) since 1D models do not accurately predict regional phase travel times and their use with regional data results in degraded event locations when combined with teleseismic phases (Bondár et al. 2004b; Myers et al. 2010).

The RSTT (Regional Seismic Travel Time) method is a seismic velocity model and computer software package that captures the major effects of three-dimensional (3D) crust and upper mantle structure on regional seismic travel times (Myers et al. 2010). RSTT was designed to be incorporated into real-time event location systems where travel times must be calculated in milliseconds on conventional computer hardware and to work seamlessly with travel time predictions for teleseismic P-waves that are based on standard 1D models (e.g., iasp91 (Kennett and Engdahl 1991), ak135 (Kennett et al. 1995)). The largest improvement in location accuracy is achieved when RSTT is used with a regional network, but improved accuracy is also achieved when RSTT is used with global networks like the IMS where at least a few stations are likely to be within regional distance of any given event. RSTT is particularly useful for low-magnitude events where signals are more reliably detected at regional stations and those stations tend to close azimuthal gaps in network coverage, thus, enabling regional stations to complement a teleseismic monitoring network. Here we describe updates to the RSTT model using a refined data set of Pn, Pg, Sn, and Lg phases.

2 The Regional Seismic Travel Time Tomography Model

The RSTT model and software are developed by the three United States National Nuclear Security Administration (NNSA) National Laboratories (Los Alamos National Laboratory, Lawrence Livermore National Laboratory, Sandia National Laboratories) in order to more accurately predict travel times from regional seismic phases (≲ 15–18°) that typically cause degradation in event location accuracy when combined with teleseismic phases (≳20°).

The RSTT model is parameterized as a global tessellation of nodes with P- and S-phase velocity profiles at each node (Fig. 1). The velocity profiles are interpolated at each node to render a 3D crust that overlies a laterally variable representation of the upper mantle velocity. The mantle velocity parameterization is specified by the velocity at the crust-mantle boundary and a linear velocity gradient as a function of depth. The crust is defined as a series of seven constant-velocity layers at each node profile: water, sediment1, sediment2, sediment3, upper crust, middle crust, and lower crust. The model also includes an additional layer that is used for calculating Pg and Lg travel times. This additional layer represents bulk effective wave speed for these phases which propagate horizontally through the crust. For a full discussion of the RSTT methodology and tomography of the Pn phase, see Myers et al. (2010).

Fig. 1
figure 1

Example of the RSTT global tessellation at ~ 1°. Inset shows more detail of tessellation with one point highlighted to show the P-wave velocity profile with depth. Dashed and solid profile lines indicate constant mantle velocity, and a gradient with depth, respectively

The global tessellation uses the GeoTess framework (Ballard et al. 2016b), an open-source code developed by Sandia National Laboratories to store spatial information in both geographic and radial dimensions. The framework is comprised of two-dimensional (2D) triangular tessellations of a unit sphere with 1D radial arrays of nodes associated with each vertex of the 2D tessellations [https://www.sandia.gov/geotess].

2.1 Updates to the RSTT Tomography Model

The original RSTT model was defined with the CRUST2.0 model (Bassin et al. 2000) for the crust overlying the upper mantle velocity and gradient. Since last publication of the RSTT methodology, the crust has been updated to the CRUST1.0 model (Laske et al. 2013) which allows more accuracy in the crust compared to the previous crustal model. In addition, for the rstt201404um version, the crustal model for the South American continent has been further updated based on a study by Assumpção et al. (2013) (hereafter, the full modified CRUST1.0 model is referred to as CRUST1.0a).

The tomography formulation outlined in Myers et al. (2010) specifically focuses on the Pn phase and its travel time calculations for the crustal and mantle portions of the assumed ray path. This formulation is also consistent for the Sn phase. The calculation of the Pn/Sn travel time is based on the Zhao and Xie (1993) method with a spatially varying mantle velocity gradient (Phillips et al. 2007), along with crustal legs on the source (α) and receiver (β) portions of the assumed ray path, and a correction for the upper mantle gradient (γ) which is always a reduction (i.e., negative value) to the travel time from the upper mantle head-wave calculation (Eq. 1) (Fig. 2) (Myers et al. 2010). The full Pn/Sn travel time (TT) is defined as

$$ TT = \sum\limits_{i = 1}^{N} {d_{i} s_{i} + \alpha + \beta + \gamma } $$
(1)
Fig. 2
figure 2

Modified from Myers et al. (2010). Xm is the length of head-wave portion of the ray path (km). di is the segment length (km) of the head-wave for calculating the total travel time. α is the event side of the path, β is the station side of the path

Model cross section and representation of ray path for Pn and Sn phases.

where d and s are the distance and slowness (i.e., inverse upper mantle velocity) in each of the i segments of the head-wave portion of the ray path, α and β are defined above, and γ is defined in Myers et al. (2010) and is related to the upper-mantle gradient.

Description of the method for computing travel times for Pg and Lg phases comes from Myers et al. (2007a). The RSTT calculation of the crustal Pg and Lg ray paths and travel times is similar to the classic definition of a P* or S* phase (Storchak et al. 2003), where the ray path is from the event to a mid-crustal layer, horizontally along the mid-crustal layer, then up to the station. The difference is that we insert a mid-crustal layer with the sole purpose of characterizing the average P-wave and S-wave velocity of horizontally propagating Pg and Lg phases, respectively. This approach allows us to accurately compute travel times for phases that actually propagate throughout the crustal waveguide and are difficult to characterize using a single ray. The Pg and Lg mid-crustal layer is not used in the calculation of Pn or Sn travel times that propagate more vertically through the crust and can be characterized by a single ray. Our formulation includes a source and receiver crustal leg component of travel time (Eq. 2) (Fig. 3). The total travel time for a crustal phase (T) is defined as

$$ T = \sum\limits_{i = 1}^{N} {x_{i} s_{i} } + \alpha^{c} + \beta^{c} $$
(2)
Fig. 3
figure 3

source and receiver station sides of the ray path, respectively

Model cross section and representative ray path for Pg and Lg phases. αc and βc represent the

where x and s are the distance and slowness (i.e., inverse of the mid-crustal velocity) in each of the i segments of the crustal head-wave portion of the ray path, and αc and βc are the source and receiver crustal leg portions of the travel time, respectively.

The separate crustal layer for Pg and Lg is only used for the calculation of those phases and not used for the calculation of Pn or Sn crustal leg components of travel time. This distinct structure for the crustal phases is used as ray paths for Pg and Lg are not well defined in ray theory and the authors chose not to mix tomographic definitions between the crustal and upper mantle phases.

2.1.1 Tomography and Validation Data Set

The data set used for tomography begins with input of various seismic bulletins compiled by the above-listed national laboratories and arrival time picks made by researchers at the national laboratories. The primary sources of bulletin data are the International Seismological Centre (ISC), including the IASPEI Reference GT List (Bondár et al. 2004a, 2008; Bondár and McLaughlin 2009; International Seismological Centre 2019); the International Data Centre (IDC) at the Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO PrepCom); the U. S. Geological Survey (USGS) National Earthquake Information Center; and several regional and local bulletins across the globe. Additional data contributions come from participants of five RSTT regional workshops sponsored by the CTBTO PrepCom, as well as temporary field campaigns and published studies, e.g. Array Network Facility (i.e., USArray) deployments in North America [https://www.usarray.org/public], the Deep Seismic Sounding (DSS) lines in Russia (Egorkin and Pavlenkova 1981; Kozlovsky 1990) that are associated with the Peaceful Nuclear Explosion (PNE) series (Table 1). Data from each of these sources are merged into a single database through a process of reconciling common events, stations, and picks (Flanagan et al. 2009). After bulletin reconciliation and preliminary event relocation, the accuracy of event hypocenter parameters is assessed for event locations based on seismic data using network coverage criteria (e.g., Bondár et al. 2004a) and published studies that determine event locations based on first-person accounts of explosives detonation (e.g., Kozlovsky 1990). A subset of events with the most associated arrival times and most accurately determined hypocenter parameters is identified for further processing.

Table 1 Seismic arrival data sources

Data processing follows the procedures outlined by Myers et al. (2011) and Simmons et al. (2012), which utilizes the Bayesloc multiple-event relative relocation algorithm (Myers et al. 2007b, 2009). This algorithm uses a non-linear Markov-Chain Monte Carlo method to simultaneously relocate events, determine the probability that phase labels are properly assigned, and estimate travel time corrections and data uncertainties. In order to minimize bias, we apply prior constraints to hypocenter parameters based on assessments of accuracy determined in the bulletin reconciliation process. Arrival times with posteriori phase-label probability greater than 0.98 are retained. Retention levels for Pn, Pg, Sn, and Lg are 90%, 93%, 73%, and 61%, respectively. Posteriori assessments of hypocenter uncertainty and data uncertainty are combined to form a datum-specific uncertainty for each tomographic ray. After Bayesloc processing the tomographic data set includes 112,168 distinct events (Fig. 4) and 2.9 million regional-phase arrivals, with a range of reported body-wave magnitudes from 0.7 to 7.5 with an average of 4.8 ± 0.7.

Fig. 4
figure 4

(top) Subset of events (white) from the full tomographic data set in which there was a minimum of three regional phases (i.e., Pn, Pg, Sn, Lg) from the IDC Reviewed Event Bulletin (REB). This subset had preferential weighting over the remaining data during tomographic inversion. Colored circles (indicating count of regional phases) represent validation events (30% of data set), randomly selected from the REB subset to leave out of the tomographic inversion and use for location testing. Validation events are sorted so events with fewer available regional phases are plotted on top. There are hundreds of events in Europe with greater than ten available phases. (bottom) Distribution of reported body-wave magnitudes (68,415) for tomography data set

For all regional phases, the tomographic starting model is based on the CRUST1.0a model, with the crustal values for Pg and Lg layers transferred from the 2014 version of the RSTT model (rstt201404um) that is publicly available on the RSTT web site [https://www.sandia.gov/rstt]. The rstt201404 version of the model uses CRUST2.0 but the CRUST1.0a model in the rstt201404um version does not include the specially-defined Pg and Lg layer values. The RSTT model framework has special Pg and Lg layers used for calculation of those phase travel times which is not consistent with CRUST1.0a, so they had to be added. The starting values for the upper mantle gradient are set to constant values of 0.0018 s−1 for both Pn and Sn and is based on past experience. There is also slight trade-off with upper mantle velocity, depending on path length, and we want to allow the data to control the upper mantle gradient. The previous rstt201404um model has a starting value of 0.0005 s−1 in areas of zero data coverage with the same 0.002 s−1 value above in areas with coverage. The average gradient for rstt201404um in areas of coverage is ~ 0.002 s−1, hence the decision to set that as the default for the current model.

The tomographic inversion uses the LSQR inversion algorithm (Paige and Saunders 1982), including both damping and smoothing (regularization) parameters that are optimized for the data set. The input event locations and origin times are not adjusted during tomography. Deeper mantle events are permitted in subduction zone areas (610 events have depths > 50 km, but the majority of events have depths < 50 km (Fig. 5)). Before the Bayesloc relocation process, events in the mantle were removed to ensure the validity of Pg and Lg phases and to better adhere to the upper mantle gradient and turning-point depth assumptions (Fig. 2). Depths were allowed to vary in the Bayesloc relocation procedure, resulting in a small percentage of events in the upper mantle. We note that the RSTT algorithm is defined for events below the Moho, but only for Pn and Sn phases (Myers et al. 2010). The tomographic inversion process solves for changes to the upper mantle slowness, the upper mantle gradient, the middle crustal values for Pg and Lg phases, and multipliers to the full crustal stack for P or S phases. The crustal stack modifier is solved for instead of standard station terms to fully incorporate slight changes to the crust into the model without the need for an external station term file, thereby permitting easier use of the model in practice. The inversion process is performed iteratively until the change in root-mean-square value (RMS) of the travel time residuals becomes small (< 1%). The inversions for various parameters are done individually at first, with damping parameters set higher than normal to minimize changes. Each parameter is held fixed while the next parameter is solved, in the order: upper mantle slowness, upper mantle gradient, crustal slowness, and crustal stack multiplier. In the final step, all parameters are solved for simultaneously, with optimized damping and smoothing values, allowing the inversion to balance the parameter fits. P- and S-phase inversions are performed separately, with the P-phase performed first and held fixed while the S-phase inversion is calculated in a similar manner to the P-phase inversions.

Fig. 5
figure 5

Histogram of event depths in the RSTT data set. The black bar on the right shows the number of events with depth > 50 km (610). Deeper mantle events were permitted in the tomography for subduction zones areas

The initial goal for an updated RSTT model is to optimize the results for use with the IMS stations. To that end, we identify phase arrivals specifically from the IDC Reviewed Event Bulletin (REB) that are associated with the tomographic data set events above to use for preferential weighting during the tomography procedures, as well as for later validation steps. Of the 112,168 tomographic data set events above, 7638 of them have at least three regional phases listed in the REB. To allow for validation events after tomographic inversion, we randomly select 30% of the events to leave out of the tomography data set, resulting in a sufficient number of validation events and regional phases (2309) (Fig. 4).

2.1.2 Pn/Pg Tomography Results

The Pn phase tomography results are the most robust, with ~ 1.9 million ray paths for the full tomographic inversion (Fig. 6). This new Pn data set is significantly increased from the ~ 600,000 Pn rays from the original Eurasian model described in Myers et al. (2010). North America, Eurasia, and Japan have the greatest ray path density, reaching over 105 ray paths in some areas for the full inversion. The weighted RMS error (by the arrival pick error) progresses from 1.57 s using the starting model to 1.34 s after three iterations of tomographic inversion. The relatively small improvement of the RMS value in the full inversion is a result of using the last RSTT model for gradient, the use of the more accurate CRUST1.0a for upper mantle velocity as the starting model, and the fact that Pn is the predominant phase in the RSTT model results. The results of the Pn tomography for upper mantle velocity and gradient are shown in Fig. 7 along with views from both the previous model (rstt201404um) and the starting model for comparison. Slower velocities are clearly associated with known areas of high seismicity/tectonic complexity, such as the Juan de Fuca Ridge area off the west coast of North America, the western United States, the mid-Atlantic Ridge/Iceland in the north Atlantic, parts of Europe, and especially along the subduction zones of Japan and southeast Asia. Cratonic areas have generally faster upper mantle velocities as expected. Upper mantle gradient values generally reflect the starting value of 0.0018 s−1, with significant departures, both low and high, near tectonically-complex areas such as the subduction zone in Japan.

Fig. 6
figure 6

Pn ray path density (log10 of the hit count) for the full RSTT tomographic inversion. Density is greatest in North America, Eurasia, and Japan, with the lowest density in the southern hemisphere

Fig. 7
figure 7

Pn tomography results for upper mantle velocity (km/s) (bottom-left) and upper mantle gradient (s−1) (bottom-right) for the RSTT tomographic inversion. For comparison, velocity for the CRUST1.0a starting model is shown at middle-left (note the presence of low-velocity mid-ocean ridge areas) with the previous velocity and gradient model (rstt201404um) at the top-left and top-right, respectively. The starting gradient value for the current tomographic inversion was 0.0018 s−1. Areas with no tomographic resolution are dimmed and are unaltered from the starting model. Note the slower velocities along mid-ocean ridges come from the values in the upper mantle for the CRUST1.0a model (Laske et al. 2013). Perturbations from the starting model are small, with more noticeable changes in the western United States and eastern Asia. Changes in gradient values are different due both from the use of the Bayesloc data as well as different starting values for the background

Ray path coverage for Pg (Fig. 8) (~ 529,000 rays) indicates much more limited coverage than for the Pn phase. Coverage is concentrated in the northern hemisphere, particularly in North America and Eurasia. There are only isolated spots of coverage in the southern hemisphere. Tomography RMS values for Pg (Fig. 9) decrease from 1.30 to 0.78 s (40% reduction) in the five iterations (after the three iterations above for Pn tomography to converge, Pn tomography was turned off). Tomography results show velocity patterns suggesting a strong correlation with known tectonic activity, such as western North America, Yellowstone National Park, Alaskan subduction zones, Southern Europe near Italy and Greece, the Baltic countries, parts of the Iranian Plateau, and parts of Japan. All of these complex areas show relatively fast velocity values compared to surrounding regions. The Pg crustal velocities are much slower however, in the complex area of Tibet where the Indian continent is colliding with Asia.

Fig. 8
figure 8

Ray coverage for Pg (log10 hit count) for the RSTT tomographic inversion. Pg coverage is far more limited than Pn, with concentrations in North America and Eurasia and very isolated spots of coverage in the southern hemisphere

Fig. 9
figure 9

Pg tomography results for the separate mid-crustal layer for the RSTT tomographic inversion (right), compared to the starting/previous model (rstt201404um, left). Areas with no tomographic resolution are dimmed and are unaltered from the starting model. Ray path coverage is minimal in the southern hemisphere, with small areas of sufficient ray path coverage for tomography. Velocity patterns suggest correlation with known tectonic areas

The crustal stack slowness multiplier result for tomography using both Pn and Pg phases indicates minimal change to the average vertical slowness of the crustal stack globally (Fig. 10). The multiplier in California and Nevada is greater than 1, indicating a slowing of the crustal stack values, while the Great Plains region just east of Yellowstone suggests faster crustal stack velocities. In Eurasia, crustal stack velocities are slightly faster for Greece, western Turkey, the Kashmir-Karakoram area and the Qaidam Basin.

Fig. 10
figure 10

Crust stack slowness modifier (or multiplier) for Pn and Pg phases for the RSTT tomographic inversion. Minimal changes are observed after tomography. Areas with no tomographic resolution are dimmed and are unaltered from the starting model (i.e., multiplier remains at 1)

2.1.3 Sn/Lg Tomography Results

Sn ray coverage (Fig. 11) shows similar patterns to the Pn ray coverage, with ~ 334,000 rays available. Coverage is concentrated in the northern hemisphere, with some of DSS lines being more obvious in the ray coverage pattern (than for Pn). Tomography RMS residuals reduce from 4.20 s with the starting model to 2.67 s in the final, a 36% decrease after five iterations. The final upper mantle velocity and gradient values are shown in Fig. 12. Upper mantle velocity patterns appear more complex than for Pn, but still tend to be associated with complex tectonic areas. Of note are the slower values in the western United States and faster values in general in the eastern United States, slower values near the Red Sea and Afghanistan, and slow values in the Yellow Sea and Japan. Upper mantle gradient values for Sn appear relatively muted around the globe, with the tomographic inversion resulting in almost no gradient in order to fit the majority of Sn travel times. This pattern is also consistent with differences between Pn and Sn amplitude variability (Phillips et al. 2016) and surface wave/receiver function studies (Nyblade et al. 2012). Two areas of relatively high Sn gradient are the western United States and Japan.

Fig. 11
figure 11

Ray coverage for Sn (log10 hit count) for the RSTT tomographic inversion. We observe a similar overall coverage pattern to Pn

Fig. 12
figure 12

Sn tomography results for upper mantle velocity (km/s) (bottom-left) and upper mantle gradient (s−1) (bottom-right) for the RSTT tomographic inversion. For comparison, velocity for the CRUST1.0a starting model is shown at middle-left (note the presence of low-velocity mid-ocean ridge areas) with the previous velocity and gradient model (rstt201404um) at the top-left and top-right, respectively. The starting gradient value for the current tomographic inversion was 0.0018 s−1. Areas with no tomographic resolution are dimmed and are unaltered from the starting model. Velocity patterns are more complex than for Pn, but appear to still be associated to tectonically complex areas. Sn gradient values are mostly minimal and consistent for large areas. Changes in gradient values are different due both from the use of the Bayesloc data as well as different starting values

Ray coverage for the Lg phase is shown in Fig. 13. There are ~ 131,000 Lg rays concentrated in the northern hemisphere, especially the United States and the Iberian Peninsula. To keep the Lg tomography independent of the Pg results, we set the starting Lg model to a constant 3.6 km/s, that is used in the iasp91 (Kennett and Engdahl 1991) and ak135 (Kennett et al. 1995) 1D velocity models. Lg tomography RMS values reduce from 4.43 to 3.81 s, a 14% reduction in RMS after seven iterations (like for Pg/Pn, Sn tomography was stopped after the five iterations above). Lg tomography results (Fig. 14) show mixed fast/slow values in the United States, with slightly more intense velocity changes in Europe and into the Middle East, such as the Iranian Plateau. The crustal stack multiplier for combined Sn/Lg phases (Fig. 15) indicates slight changes from the default of 1 (i.e., no change), but with higher amplitude values than those for the Pn/Pg tomography. The values for Italy suggest slower crustal stack values, with the surrounding mountains suggesting faster crustal stack values.

Fig. 13
figure 13

Lg ray coverage (log10 hit count) for the RSTT tomographic inversion. Like Pg to Pn, Lg is more limited in coverage than Sn but has significantly denser coverage than Pg, particularly in the northern hemisphere

Fig. 14
figure 14

Lg tomography results for the RSTT tomographic inversion (right), compared to the previous model (rstt201404um, left). The starting model for Lg was set to 3.6 km/s to assure independence from Pg tomography results. Areas with no tomographic resolution are dimmed and are unaltered from the starting model

Fig. 15
figure 15

Crustal modifier (or multiplier) values for combined Sn and Lg phases for the RSTT tomographic inversion. Areas with no tomographic resolution are dimmed and are unaltered from the starting model. The more significant changes to the default of 1 are observed in Europe and Japan

3 Validation of RSTT Tomography Model

Validation events are selected from the full Bayesloc data set where the events are listed in the IDC REB and have a minimum of three regional phases (i.e., Pn, Pg, Sn, Lg) listed in the REB (Fig. 4). The event hypocenters are set to be those determined during the Bayesloc relative relocation procedure, with validation events limited to those with depths < 50 km. All validation events are assumed to be accurate to 5 km or better (i.e., GT5) based on uncertainty estimates after Bayesloc relocation. All regional phases are set to “defining” whether or not they are originally defining in the REB. The arrival times are those directly reported in the REB.

For relocation of the 2309 validation events, we use the LocOO3D algorithm developed by Sandia National Laboratories [https://www.sandia.gov/salsa3d/Software.html] (Ballard et al. 2008, 2009) which permits using standard 1D travel-time tables (e.g., iasp91, ak135), source-specific correction surfaces (SSSCs) for each IMS station and phase (Firbas et al. 1998) (converted to GeoTess format), ray-bending through 3D models, and full RSTT models for travel-time prediction. All depths are fixed at the Bayesloc result to reduce the trade-off between depth and origin time. We solve for a 90% coverage ellipse (assuming a priori estimates of model error—1D or 2D), using only regional REB phases (Pn, Pg, Sn, Lg), with a limit of 17 degrees to be consistent with the RSTT assumptions for turning point depths above the 410 discontinuity. No phases are removed during relocations because of large residuals, and all primary and auxiliary IMS stations are used that are reported in the REB.

We compare epicenter mislocations after phase-specific and full regional phase single-event relocations to the Bayesloc-defined locations. Mislocation results are binned with a one-arrival overlap. For velocity models, we use (1) the 1D iasp91 model (the 1D model used by the IDC), (2) SSSCs produced by the IDC (IDC-SSSCs), and (3) the RSTT model developed above (with validation events removed before tomographic inversion) using a standard 1D, distance-dependent model uncertainty. For the IDC-SSSC relocations, only about one-third of the primary/auxiliary stations have defined SSSCs (Fig. 16); the remaining stations will default to the iasp91 model. Due to this limitation of existing IMS stations with defined IDC-SSSCs, we compare mislocations when using the full IMS network (primary + auxiliary) and when using only those IMS stations with defined IDC-SSSCs. The companion paper, “Updates to the Regional Seismic Travel Time (RSTT) Model: 2. Path-dependent Travel-time Uncertainty”, this issue, describes validation results for coverage ellipse size using the models listed above and including path-dependent travel-time uncertainties for the RSTT model.

Fig. 16
figure 16

source-specific stations corrections (IDC-SSSCs) defined for regional phases (Pn, Pg, Sn, Lg) (Firbas et al. 1998)

Map of International Monitoring System primary (red) and auxiliary (blue) seismic stations. Circles indicate which stations had original IDC

3.1 Mislocation Results

Figure 17 shows validation events when selecting only specific regional phases. Of the 2309 possible REB validation events, 1395 had a minimum of three available arrivals for Pn-only relocations for all IMS stations and 914 when considering only stations with defined IDC-SSSCs. Figure 18 shows the single-event Pn mislocation values for relocations with the velocity models specified above. The RSTT model shows significant reduction in mislocation (14.7 km) when compared to both the iasp91 model and the IDC-SSSCs (18.6 km and 18.2 km, respectively). When considering only stations with defined IDC-SSSCs, the median mislocation values for RSTT are very similar to those when using all IMS stations with the median value increasing to 18.7 km due to the distribution relative to the number of arrivals. SSSC-only relocations for both iasp91 and IDC-SSSCs show an increase in median mislocation compared to their values when all IMS stations are used until about seven total arrivals when the values intersect those for all stations.

Fig. 17
figure 17

Validation events with a minimum of three phase-specific REB arrivals compared to all possible REB events (from tomography data set) with a minimum of three phase-specific REB arrivals (top-left) 1395 Pn-only validation events (blue) with 4533 Pn-only tomography events (white); (top-right) 62 Pg-only validation events (cyan) with 169 Pg-only tomography events (white); (bottom-left) 528 Sn-only validation events (red) with 1612 Sn-only tomography events (white); (bottom-right) 434 Lg-only validation events (green) with 1409 Lg-only tomography events (white)

Fig. 18
figure 18

Validation results for Pn-only relocations using the velocity models specified and all IMS stations (solid lines) or only those stations having defined IDC-SSSCs (dashed lines). (top) Median mislocation relative to the number of Pn arrivals. Overall median mislocation for each velocity model is specified in the legend. The RSTT model (blue) shows significant reduction in median mislocation values. (bottom) Histogram of the number of events with specific arrival counts for the iasp91 model using all IMS stations (black) or using only those stations with defined IDC-SSSCs (gray) to demonstrate numbers of arrivals in each bin

Of the 2309 possible REB validation events, only 62 had a minimum of three available arrivals for Pg-only relocations when using all IMS stations and 39 when using only stations with defined IDC-SSSCs. The number of available Pg arrivals is the lowest of any of those available for the other regional phases (i.e., Pn, Sn, Lg). Figure 19 shows the single-event relocation results. The results when using the RSTT model are higher (17.5 km) compared to the iasp91 model, having the lowest median mislocation of 15.8 km. The IDC-SSSCs have significantly higher median mislocation values for the Pg-only events (44.1 km). When restricting to only stations with defined IDC-SSSCs, results using iasp91 do not change significantly, whereas RSTT results do degrade slightly to 22.7 km and IDC-SSSC results flatten with the number of arrivals with a median mislocation value increasing to 46.8 km. The reason that using the iasp91 travel times results in the most accurate Pg-based locations may be due to picking procedures for this commonly emergent phase that are guided by predicted phase arrival times, but this result deserves further study.

Fig. 19
figure 19

Validation results for Pg-only relocations using the velocity models specified. Features of plots are similar to Fig. 18 above. The mislocation result for the RSTT model (17.5 km) is greater than for the iasp91 model (15.8 km), with both models having significantly lower mislocations compare to the IDC-SSSCs (44.1 km). When restricting to only stations with original IDC-SSSCs, the IDC-SSSC mislocation line flattens but results in increased median mislocation of 46.8 km. The RSTT results when using only stations with IDC-SSSCs degrade to a median of 22.7 km

There are 528 REB validation events that had three available arrivals for Sn-only relocations when using all IMS stations and 255 when restricting to only stations with defined IDC-SSSCs. Figure 20 shows the single-event relocation results for Sn-only validation events. As for the Pn-only relocations (Fig. 18), the median mislocation results again show that the RSTT model significantly reduces mislocation (19.9 km) over the variable number of phases compared to the iasp91 model (26.2 km) or the IDC-SSSCs (25.5 km). When considering relocations using only stations with defined IDC-SSSCs, median mislocations for iasp91, IDC-SSSC, and RSTT models are 34.3, 27.2, and 21.9 km, respectively. The relative differences are similar between the models at 3–4 available Sn arrivals. From 5–7 Sn arrivals, the IDC-SSSC and RSTT models have similar median mislocation values.

Fig. 20
figure 20

Validation results for Sn-only relocations using the velocity models specified. Features of plots are similar to Fig. 18 above. The RSTT model shows reduced mislocation values (19.9 km) compared to the iasp91 (26.2 km) and IDC-SSSCs (25.5 km). Results when restricted to only stations with defined IDC-SSSCs show similar relative pattern for 3–4 available Sn arrivals. From 5 to 7 Sn arrivals, the IDC-SSSC and RSTT models have similar results

There are 434 REB validation events that had three available arrivals for Lg-only relocations and 319 events when stations are restricted to those with defined IDC-SSSCs. Figure 21 shows the single-event relocation results for Lg-only validation events. The RSTT model shows a significant reduction in median mislocation (16.2 km), compared to other tested models. The results for iasp91 and the IDC-SSSCs all show similar median mislocation values at 18.9 and 18.7 km, respectively. There are significant fluctuations in median values with increasing number of arrivals for all the model versions tested. Since Lg arrival picks tend to be the most variable of all the regional phases, this result is not surprising. When restricting to stations with defined IDC-SSSCs, median mislocations values are similar for iasp91 and RSTT from 4 to 5 available Lg phases, with the IDC-SSSC results being relatively degraded, and diverge more from 6 to 7 phases, with all models showing an improvement compared to using all IMS stations.

Fig. 21
figure 21

Validation results for Lg-only relocations using the velocity models specified. Features of plots are similar to Fig. 18 above. For the RSTT and IDC-SSSCs, median mislocation values vary considerably over the number of arrivals. The median RSTT values are the lowest of all the models compared

The above relocation tests use only one phase at a time, attempting to compare how the RSTT model for each phase would affect seismic locations. In order to test “normal” regional events, we allowed any combination of regional phases to be used. Figure 22 shows the percent of each phase for each REB validation event. The Pg phase is usually no higher than 30% of the total phases for an event, with the Pn phase being the most predominant. The Sn and Lg phases can get as high as ~ 50%.

Fig. 22
figure 22

Percent of each regional phase per event for testing relocation results with all available regional phases. Gray circles indicate there are no phases of the indicated name for that event

Figure 23 shows the relocation metric results when using any available REB regional phase. The number of the possible 2309 validation events with a minimum of any three regional phases drops from 1779 when using all IMS stations to 1149 when restricting to stations with defined IDC-SSSCs. The median mislocation results show the overall reduction in mislocation when using the RSTT model is 15.3 km. The IDC-SSSCs reach 18.9 km, while the iasp91 model has a median value of 22.1 km. When considering only the results using stations with defined IDC-SSSCs, all models show a decrease in median mislocations for ~ 3–6 arrivals compared to using all IMS stations and then an increase in median mislocation from ~ 6 to 10 arrivals. Beyond 10 arrivals the models have similar results with RSTT displaying slightly lower mislocations overall in this arrival range. We compare the overall bias (mislocation in a particular direction) of each model by plotting mislocation vectors for each validation event (Fig. 24). Each model results in some extreme mislocation values (> 60 km, mostly related to a small number of available phases (i.e., 3–4)) that are readily apparent, the use of the RSTT model does show a reduction in the number of extreme mislocations (209) over the iasp91 (308) and IDC-SSSC (256) models as well as more areas of reduced overall model bias. Particular areas where this is the most apparent include Alaska, where the iasp91 and IDC-SSSCs produce more elongated mislocation vectors in the northeast-southwest directions, the western United States, the north Atlantic, central Asia, the eastern Asian subduction zones, and the south Pacific regions. The IDC-SSSCs do show a clear reduction in mislocation for a cluster of events in the central United States.

Fig. 23
figure 23

Validation results for All Regional Phase relocations using the velocity models specified. Features of plots are similar to Fig. 18 above. The RSTT model shows significantly lower median mislocation values (15.3 km) compared to iasp91 (22.1 km) and IDC-SSSCs (18.9 km)

Fig. 24
figure 24

Mislocation vectors for validation event relocations using all regional phases and all IMS stations with models (1) iasp91 (top), (2) IDC-SSSCs (middle), and (3) RSTT (bottom). Events shown are GT locations from Bayesloc relocation process with colors indicating the number of phases used for relocation, vectors indicating direction of the relocated epicenter using the model listed, and vector length for amount of mislocation (exaggerated, see key). Overall the RSTT model results in fewer extreme mislocations (> 60 km) as well as fewer instances of localized bias. The IDC-SSSCs do reduce bias for a localized set of events in the central United States compared to iasp91 and RSTT. Relocations for iasp91 and IDC-SSSCs will be identical for events in the southern hemisphere due to the lack of stations with IDC-SSSCs

4 Discussion and Conclusions

The RSTT tomographic method and results are carefully determined in order to allow for the assumptions concerning ray turning-point depths. The calculated turning-point depth follows from the assumption of a starting upper mantle velocity and constant gradient, with any ray being rejected showing a significant deviation from this assumption (i.e., turning point depth combined with upper mantle gradient value being much less than a threshold value). The development of the Bayesloc data set, which allows a consistent quality control of the event locations as well as the data, also allows for further checks on the validity of the RSTT mantle velocity-gradient assumption and for fewer data outliers. The RSTT software itself does check for violations of the mantle velocity-gradient assumption as well as violations of ray assumptions for short ray paths (≲ 1 deg) and inconsistencies with any resulting travel times calculated from the model during tomography.

The resulting tomography model has all four regional phases included and has ray coverage that allows for velocity estimates for a large part of the Earth, including near mid-ocean ridges, as well as island stations. The use of the CRUST1.0a model as the starting model for tomography is clearly evident in the continued patterns of slower velocities along the mid-ocean ridges, in the western United States, Europe, and the east Asian subduction areas. The small changes in upper mantle velocity values from the starting to final Pn tomography model suggests that the CRUST1.0a model is a fairly accurate Pn velocity model for the upper mantle. The changes in the Sn velocity values from the starting CRUST1.0a to the final model are more pronounced than for Pn, particularly indicated by the increase in higher velocities in Europe into western Asia. Overall, the pattern of Pn and Sn velocities is similar between the last RSTT model version (rstt201404um) and the current tomography inversion values, barring the lack of mid-ocean ridge patterns in the previous model, although there is a suggestion of those features in the previous model where there was sufficient data coverage (similar to current model). Velocity and gradient values are consistent with global/regional tectonic structures, with the highest gradient values being associated mostly with subduction zones. Gradient values do vary significantly between the current and previous RSTT models, as suggested by values in central and southeast Asia. The current tomographic inversion employs slightly heavier constraints on gradient values (i.e., damping and smoothing). Crustal velocities (i.e., Pg, Lg) are also consistent with tectonic areas, although some faster velocities occur near the subduction zones and it is unclear whether these are from actual crustal variations or due to mislabeled mantle phases. Efforts are made prior to tomography to remove data exhibiting more obvious phase naming issues. The Pg velocity values tend to be fairly consistent with the previous model version. The Lg inversion has notable differences compared to the previous/starting model for the United States and for most of Europe and western Asia. When the rstt201404um model version was produced, it did not include the full US Array data set coverage that is included in the current model. The current Lg inversion did not use the Pg velocities with a Vp/Vs ratio for the starting model as was done for the rstt201404um version. This kept the Lg results independent of Pg, resulting in more obvious Lg velocity differences than the previous version. This is readily apparent in Europe (e.g. eastern Hellenic Arc), where velocity differences from rstt201404um tend to correspond to areas with a paucity of Lg ray paths (Fig. 13).

Validation efforts for the updated RSTT model attempt to compare both differences in velocity models as well as changes in the network to enable a direct comparison of the original IDC-SSSC surfaces to the new RSTT model. The Pn and Sn mislocation results show a significant reduction in median mislocation when using the RSTT model over the iasp91 model, and to a slightly lesser extent, over the IDC-SSSCs. This reduction in mislocation is consistent from the minimum number of available phases (3) to ~ 10 phases, just below the maximum observed of 12. These results are also consistent when the stations used are restricted to those with defined IDC-SSSCs. The stations that have defined IDC-SSSCs (Fig. 16) are exclusively in the northern hemisphere. The use of the IDC-SSSCs for events in the northern hemisphere (using all regional phases) shows a reduced median mislocation (17.8 km) over using iasp91 (20.6 km), while RSTT exhibits a further reduction in mislocation for northern hemisphere events (14.5 km) due to station coverage and data density available for tomography. For events in the southern hemisphere, RSTT provides a significant benefit over the current use of IDC-SSSCs for relocation, as the IDC-SSSCs result in no differences from those just using the iasp91 model. The use of RSTT with all regional phases for events in the southern hemisphere results in a significant reduction in median mislocation (42.0 km for iasp91, 30.2 km for RSTT).

The RSTT results for Pg phases do not show the significant improvement for relocation as for Pn phase validation. This could be due to a very limited number of validation events for Pg (62) and/or that the event locations mostly fall in one general area (North America). It is not unexpected for results using only IDC-SSSC stations to be slightly worse than when using all IMS stations given reduced station coverage as well as more events with only three stations used for location as compared to higher numbers of stations. The Pg validation results using the IDC-SSSC correction surfaces are surprisingly poor, for either all IMS stations and only those with SSSCs (44.1 and 46.8 km, respectively), and the root cause may be the arrival-time picking methodology that is guided by iasp91 arrival time predictions. Future validation of RSTT or other models should focus on increasing the availability of validation events for use with Pg phases.

The Lg mislocation results appear more erratic than for Pn or Sn (the lack of validation results for Pg make comparisons difficult). The Lg phase is known to be difficult to pick consistently due to its emergent nature, presumably contributing to the erratic location results. Beyond 7 available arrivals, the RSTT median values sharply increase from 7.8 to 14.5 km (from a count of 21 events to eight events, respectively) (Fig. 21). The number of events falls again to five for nine available arrivals, limiting robust statistical comparisons.

For results using all regional phases, all the tested models show reduced median mislocations for fewer than six arrivals when only considering stations with defined IDC-SSSCs, compared to using all IMS stations. This is immediately followed by an increase in median mislocation for the ~ 6–10 arrival range and for all tested models compared to when using all IMS stations. If this pattern were only observed for the IDC-SSSC and/or the RSTT relocations, it might be explained by the availability of data for tomography and/or validation for stations with defined IDC-SSSCs, presumably because the original IDC-SSSCs were developed at the time for those stations with significant data coverage. Given that the pattern is observed for iasp91 relocations as well, this suggests that the events with fewer than 10 phases have a combination of phases and stations that either improves location results (fewer than six arrivals) or degrades results (~ 6–10 arrivals). For events resulting in less than six arrivals for only defined IDC-SSSC stations, the majority have combinations of P and S phases (i.e., mostly Pn and Sn); for between six and ten arrivals, they consist of a large majority of Pn arrivals. This lack of S phases for the restricted set of stations during relocation appears to have resulted in increased mislocation values.

Overall, the updated RSTT tomography model permits the determination of more accurate travel times for events spanning most of the Earth, thus providing reduced mislocation over the standard iasp91 model and the original IDC-SSSCs, which are geographically limited to the northern hemisphere, and more confidence in event location results.