Abstract
A function of global monitoring of nuclear explosions is the development of Earth models for predicting seismic travel times for more accurate calculation of event locations. Most monitoring agencies rely on fast, distance-dependent one-dimensional (1D) Earth models to calculate seismic event locations quickly and in near real-time. RSTT (Regional Seismic Travel Time) is a seismic velocity model and computer software package that captures the major effects of three-dimensional crust and upper mantle structure on regional seismic travel times, while still allowing for fast prediction speed (milliseconds). We describe updates to the RSTT model using a refined data set of regional phases (i.e., Pn, Pg, Sn, Lg) using the Bayesloc relative relocation algorithm. The tomographic inversion shown here acts to refine the previous RSTT public model (rstt201404um) and displays significant features related to areas of global tectonic complexity as well as further reduction in arrival residual values. Validation of the updated RSTT model demonstrates significant reduction in median epicenter mislocation (15.3 km) using all regional phases compared to the iasp91 1D model (22.1 km) as well as to the current station correction approach used at the Comprehensive Nuclear-Test-Ban Treaty Organization International Data Centre (18.9 km).
Similar content being viewed by others
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).
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
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
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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%.
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.
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.
Availability of Data and Material
Seismic arrival measurements specifically from the US National Laboratories are not publicly available. Other seismic bulletin data we obtained from public sites, except for the International Data Centre Reviewed Event Bulletin which is only available to member states to the Comprehensive Nuclear-Test-Ban Treaty Organization.
Code Availability
The tomographic inversion code is a custom code from Los Alamos National Laboratory. Other algorithms are available are as specified in the manuscript.
References
Amaru, M. (2007). Global travel time tomography with 3-D reference models (Vol. 274). Geologica Ultraiectina: Utrecht University.
Assumpção, M., Feng, M., Tassara, A., & Juliá, J. (2013). Models of crustal thickness for South America from seismic refraction, receiver functions and surface wave tomography. Tectonophysics, 609, 82–96. https://doi.org/10.1016/j.tecto.2012.11.014
Ballard, S., Hipp, J., Begnaud, M., Young, C., Encarnacao, A., Chael, E., et al. (2016a). SALSA3D—A tomographic model of compressional wave slowness in the Earth’s mantle for improved travel time prediction and travel time prediction uncertainty. Bulletin of the Seismological Society of America. https://doi.org/10.1785/0120150271
Ballard, S., Hipp, J., Kraus, B., Encarnacao, A., & Young, C. (2016b). GeoTess: A generalized earth model software utility. Seismological Research Letters, 87(3), 719–725. https://doi.org/10.1785/0220150222
Ballard, S., Hipp, J. R., & Young, C. J. (2009). Efficient and accurate calculation of ray theory seismic travel time through variable resolution 3D earth models. Seismological Research Letters, 80(6), 989–998. https://doi.org/10.1785/gssrl.80.6.989
Ballard, S., Hipp, J., Young, C., Barker, G. T., & Chang, M. (2008). Implementation of a pseudo-bending seismic travel time calculator in a distributed parallel computing environment. In 30th Monitoring Research Review, Portsmouth, VA, 23–25 Sept. (pp. 338–346, paper 332–302).
Bassin, C., Laske, G., & Masters, G. (2000). The current limits of resolution for surface wave tomography in North America. EOS, Transactions of the American Geophysical Union, 81, F897.
Bondár, I., Bergman, E., Engdahl, E., Kohl, B., Kung, Y.-L., & McLaughlin, K. (2008). A hybrid multiple event location technique to obtain ground truth event locations. Geophysical Journal International. https://doi.org/10.1111/j.1365-246X.2008.03867.x
Bondár, I., Engdahl, E., Yang, X., Ghalib, H., Hofstetter, A., Kirichenko, V., et al. (2004). Collection of a reference event set for regional and teleseismic location calibration. Bulletin of the Seismological Society of America, 94, 1528–1545. https://doi.org/10.1785/012003128
Bondár, I., & McLaughlin, K. L. (2009). A new ground truth data set for seismic studies. Seismological Research Letters, 80(3), 465–472.
Bondár, I., Myers, S. C., Engdahl, E. R., & Bergman, E. A. (2004). Epicentre accuracy based on seismic network criteria. Geophysical Journal International, 156, 483–496.
Burdick, S., Li, C., Martynov, V., Cox, T., Eakins, J., Mulder, T., et al. (2008). Upper mantle heterogeneity beneath North America from travel time tomography with global and USArray transportable array data. Seismological Research Letters, 79(3), 384–392. https://doi.org/10.1785/gssrl.79.3.384
Chang, S.-J., van der Lee, S., Flanagan, M., Bedle, H., Marone, F., Matzel, E. M., et al. (2010). Joint inversion for three-dimensional S velocity mantle structure along the Tethyan margin. Journal of Geophysical Research. https://doi.org/10.1029/2009JB007204
Egorkin, A. V., & Pavlenkova, N. I. (1981). Studies of mantle structure in the U.S.S.R. territory on long-range seismic profiles. Physics of the Earth and Planetary Interiors, 25, 12–26.
Engdahl, E. R., Di Giacomo, D., Sakarya, B., Gkarlaouni, C. G., Harris, J., & Storchak, D. A. (2020). ISC-EHB 1964–2016, an improved data set for studies of earth structure and global seismicity. Earth and Space Science. https://doi.org/10.1029/2019EA000897
Engdahl, E. R., van der Hilst, R., & Buland, R. (1998). Global teleseismic earthquake relocation with improved travel-times and procedures for depth determination. Bulletin of the Seismological Society of America, 88(3), 722–743.
Firbas, P., Fuchs, K., & Mooney, W. D. (1998). Calibration of seismograph network may meet Test Ban Treaty’s monitoring needs. EOS, Transactions of the American Geophysical Union, 79, 413–421.
Flanagan, M. P., Dodge, D., & Myers, S. C. (2009). GT Merge Process: Version 2.0, Lawrence Livermore National Laboratory, LLNL-TR-410856, 21 pp.
International Seismological Centre. (2019). IASPEI Reference Event (GT) List. https://doi.org/10.31905/32NSJF7V.
Kennett, B. L. N., & Engdahl, E. R. (1991). Traveltimes for global earthquake location and phase identification. Geophysical Journal International, 105(2), 429–465.
Kennett, B. L. N., Engdahl, E. R., & Buland, R. (1995). Constraints on seismic velocities in the Earth from travel times. Geophysical Journal International, 122, 108–124.
Kozlovsky, Y. A. (1990). The USSR integrated program of continental crust investigations and studies of the earth’s deep smacmre under the globus project. In K. Fuchs, Y. A. Kozlovsky, I. Krivtsov, & M. D. Zoback (Eds.), Super-deep continental drilling and deep geophysical sounding (pp. 90–103). Berlin: Springer.
Laske, G., Masters, G., Ma, Z., & Pasyanos, M. (2013). Update on CRUST1.0—A 1-degree global model of Earth's crust. European Geophysical Union (Vol. 15, pp. EGU2013-2658). Vienna, Austria.
Li, C., van der Hilst, R. D., Engdahl, E. R., & Burdick, S. (2008). A new global model for P wave speed variations in Earth’s mantle. Geochemistry, Geophysics, Geosystems, 9(5), Q05018.
Myers, S. C., Ballard, S., Rowe, C. A., Wagner, G. S., Antolik, M. S., Phillips, W. S., et al. (2007a). Tomography and methos of travel-time calculation for regional seismic location. In 29th Monitoring Research Review: Ground-based Nuclear Explosion Monitoring Technologies, Denver, Colorado, September 25–27, 2007 (pp. 414–423).
Myers, S. C., Begnaud, M. L., Ballard, S., Pasyanos, M. E., Phillips, W. S., Ramirez, A. L., et al. (2010). A crust and upper mantle model of Eurasia and North Africa for Pn travel time calculation. Bulletin of the Seismological Society of America, 100(2), 640–656.
Myers, S. C., Johannesson, G., & Hanley, W. (2007). A Bayesian hierarchical method for multiple-event seismic location. Geophysical Journal International, 171, 1049–1063.
Myers, S. C., Johannesson, G., & Hanley, W. (2009). Incorporation of probabilistic seismic phase labels into a Bayesian multiple-event seismic locator. Geophysical Journal International, 177, 193–204.
Myers, S. C., Johannesson, G., & Simmons, N. A. (2011). Global-scale P wave tomography optimized for prediction of teleseismic and regional travel times for Middle East events: 1. Data set development. Journal of Geophysical Research, 116(B04304), 13. https://doi.org/10.1029/2010JB007967
Nyblade, A. A., Juliá, J., Matzel, E. M., Raveloson, A. H., & Rodgers, A. J. (2012). Developing regionalized models of lithosphere thickness and velocity structure across Eurasia and the middle east from jointly inverting P-wave and S-wave receiver functions with Rayleigh wave group and phase velocities. In 2012 Monitoring Research Review: Ground-based Nuclear Explosion Monitoring Technologies, Albuquerque, New Mexico, (pp. 66–74).
Paige, C. C., & Saunders, M. A. (1982). LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software, 8(1), 43–71.
Phillips, W. S., Begnaud, M. L., Rowe, C. A., Steck, L. K., Myers, S. C., Pasyanos, M., et al. (2007). Accounting for lateral variations of the upper mantle gradient in Pn tomography studies. Geophysical Research Letters, 34, L14312. https://doi.org/10.1029/2007GL029338
Phillips, W. S., Fisk, M. D., Stead, R. J., Begnaud, M. L., Yang, X., & Ballard, S. (2016). Three-component high frequency amplitude models for discrimination and yield estimation. Seismological Society of America (Vol. 87, pp. 498). Reno, Nevada
Simmons, N. A., Forte, A. M., Boschi, L., & Grand, S. P. (2010). GyPSuM: A joint tomographic model of mantle density and seismic wave speeds. Journal of Geophysical Research. https://doi.org/10.1029/2010JB007631
Simmons, N. A., Myers, S. C., Johannesson, G., & Matzel, E. (2012). LLNL-G3Dv3: Global P-wave tomography model for improved regional and teleseismic travel time prediction. Journal of Geophysical Research. https://doi.org/10.1029/2012JB009525
Steck, L. K., Begnaud, M. L., Phillips, W. S., & Stead, R. (2011). Tomography of crustal P and S travel times across the western United States. Journal of Geophysical Research. https://doi.org/10.1029/2011JB008260
Storchak, D. A., Schweitzer, J., & Bormann, P. (2003). The IASPEI standard seismic phase list. Seismological Research Letters, 74(6), 761–772.
Syracuse, E. M., Zhang, H., & Maceira, M. (2017). Joint inversion of seismic and gravity data for imaging seismic velocity structure of the crust and upper mantle beneath, Utah, United States. Tectonophysics, 718, 105–117. https://doi.org/10.1016/j.tecto.2017.07.005
Vaidya, S., Engdahl, R., Le Bras, R., Koch, K., & Dahlman, O. (2009). Strategic initiative in support of CTBT data processing: vDEC (virtual Data Exploitation Centre). Abstract DM-01/A presented at the CTBTO International Scientific Studies, Hofburg, Vienna, Austria, 10–12 June 2009.
Weston, J., Engdahl, E. R., Harris, J., Di Giacomo, D., & Storchak, D. A. (2018). Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bulletin of the Seismological Society of America, 88(3), 722–743.
Zhao, L.-S., & Xie, J. (1993). Lateral variations in compressional velocities beneath the Tibetan Plateau from Pn traveltime tomography. Geophysical Journal International, 115(3), 1070–1084. https://doi.org/10.1111/j.1365-246X.1993.tb01510.x
Acknowledgements
We thank the Office of Nuclear Detonation Detection (NA-222) as well as the Office of Nuclear Verification (NA-243) within the National Nuclear Security Administration for funding this effort. The authors wish to acknowledge the data contributions by the IDC, as well as the ISC, USGS, and many other individuals and organizations, without which the RSTT model and software could not have been developed and refined. Further acknowledgements extend to staff and researchers at Los Alamos, Lawrence Livermore, and Sandia National Laboratories who helped acquire data as well as provided input and tested various versions of models and software over the years of RSTT development. The authors thank Ronan Le Bras for providing relevant and insightful suggestions to this manuscript. This work was performed under the auspices of the U.S. Department of Energy by Los Alamos National Laboratory, Lawrence Livermore National Laboratory, and Sandia National Laboratories. Los Alamos National Laboratory is managed by Triad National Security, LLC, under contract no. 89233218CNA000001. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This article is Los Alamos National Laboratory document LA-UR-20-20966. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.
Funding
United States National Nuclear Security Administration, Office of Defense Nuclear Nonproliferation.
Author information
Authors and Affiliations
Contributions
All authors contributed to the study conception and design. The RSTT software and framework were developed by MB, SM, BY, JH, and WSP. The tomography algorithm was developed by MB, WSP, SM, and BY. Data collection and analysis was performed by DD, SM, and MB. The first draft of the manuscript was written by MB and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Begnaud, M.L., Myers, S.C., Young, B. et al. Updates to the Regional Seismic Travel Time (RSTT) Model: 1. Tomography. Pure Appl. Geophys. 178, 2475–2498 (2021). https://doi.org/10.1007/s00024-020-02619-5
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00024-020-02619-5