Updates to the Regional Seismic Travel Time (RSTT) Model: 2. Path-dependent Travel-time Uncertainty

Abstract

The regional seismic travel time (RSTT) model and software were developed to improve travel-time prediction accuracy by accounting for three-dimensional crust and upper mantle structure. Travel-time uncertainty estimates are used in the process of associating seismic phases to events and to accurately calculate location uncertainty bounds (i.e. event location error ellipses). We improve on the current distance-dependent uncertainty parameterization for RSTT using a random effects model to estimate slowness (inverse velocity) uncertainty as a mean squared error for each model parameter. The random effects model separates the error between observed slowness and model predicted slowness into bias and random components. The path-specific travel-time uncertainty is calculated by integrating these mean squared errors along a seismic-phase ray path. We demonstrate that event location error ellipses computed for a 90% coverage ellipse metric (used by the Comprehensive Nuclear-Test-Ban Treaty Organization International Data Centre (IDC)), and using the path-specific travel-time uncertainty approach, are more representative (median 82.5% ellipse percentage) of true location error than error ellipses computed using distance-dependent travel-time uncertainties (median 70.1%). We also demonstrate measurable improvement in location uncertainties using the RSTT method compared to the current station correction approach used at the IDC (median 74.3% coverage ellipse).

1 Introduction

For standard seismic event location in real-time, the fastest method for calculating travel times has been the use of one-dimensional (1D) Earth models that vary only with depth. In practice, these 1D models have been used to calculate travel-time lookup tables for various seismic phases that only require the event depth and source–receiver distance as input. In addition to travel time, these lookup tables also include 1D information for estimating distance-dependent travel-time uncertainty. However, 1D models do not accurately predict regional phase travel times with uncertainties, and their use with regional data results in degraded event locations when combined with teleseismic phases (Bondár et al. 2004; Myers et al. 2010).

Many authors have developed higher-dimensional (i.e., two-dimensional (2D) and three-dimensional (3D)) global and regional Earth models for the crust and upper mantle, 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), and for others to better predict seismic travel times for event location (e.g., Ballard et al. 2016a; Phillips et al. 2007; Simmons et al. 2012). Using Earth models to better predict travel times is a primary focus for organizations responsible for, as well as interested researchers involved in, nuclear explosion monitoring.

The need and desire to quickly locate all seismic events has pushed magnitude thresholds lower and lower. Online data access is more prevalent than ever with near real-time data and results from the Global Seismic Network (GSN) and the International Seismological Centre (ISC), for example, as well as the utilization of existing global, regional and local networks (e.g., collections by the United States Geological Survey (USGS) and the International Monitoring System (IMS)). As magnitude thresholds decrease, the reliance on more regional-distance seismic stations (i.e., ≲15–18°) is increasing in order to have sufficient identified seismic phases with which to locate the smaller magnitude events.

RSTT (Regional Seismic Travel Time) is a seismic velocity model and computer software package that captures the major effects of 3D crust and upper mantle structure on regional seismic travel times (Myers et al. 2010). This model and software package 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) and ak135 (Kennett et al. 1995)). Using RSTT results in more accurate predictions for regional-phase (i.e., Pn, Pg, Sn, or Lg) travel times, and combining these predictions in a global network like the IMS, where at least a few stations are likely to be within regional distance of any given event, also results in improved event location accuracy.

The use of Earth models for travel-time prediction also includes estimates of prediction uncertainty or “model” error. Traditionally, error estimates for higher-dimensional Earth models have been relegated to almost an afterthought, with a simple distance-dependent phase uncertainty estimate calculated for event location. A fundamental goal of monitoring agencies is to produce location uncertainty bounds that are representative of true error, and representative uncertainty estimates for travel-time calculations are key to achieving that goal. Higher-dimensional Earth models should include robust estimates of their prediction uncertainties in order to provide accurate size (i.e., precision) and orientations of error ellipses, especially considering that a maximum area of 1000 km² is specified in the Comprehensive Nuclear-Test-Ban Treaty (CTBT) for any possible onsite inspection [https://www.ctbto.org/fileadmin/content/treaty/treatytext.tt.html].

Here we describe the formulation of global, path-dependent, travel-time uncertainty estimates for the updated RSTT tomography model using a refined data set of Pn, Pg, Sn, and Lg phases (see companion paper “Updates to the Regional Seismic Travel Time (RSTT) Model: 1. Tomography”, this issue).

2 The Regional Seismic Travel Time Model

The RSTT model and software was 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°). For a full description of the RSTT model parameterization, framework, data set, and tomography results, see companion paper “Updates to the Regional Seismic Travel Time (RSTT) Model: 1. Tomography”, this issue.

The tomography data set is a result of relocating a reconciled global bulletin that includes the ISC, IDC Reviewed Event Bulletin (REB), regional-network, and national laboratory data. The Bayesloc relative relocation algorithm (Myers et al. 2007, 2009) was used to improve event locations and data quality and consistency, following a similar event clustering procedure described by Simmons et al. (2012). This procedure utilizes stochastic constraints on locations from previous studies, including known explosion locations and locations determined using dense local networks that meet the “ground-truth” (GT) criteria of Bondár et al. (2004), to minimize location bias and simultaneously assess data consistency and quality. The full bulletin data for input to Bayesloc includes events whose epicenter locations were determined to be GT0–5 km (i.e., between 0 and 5 km from ground truth) based on a priori knowledge (e.g., nuclear and chemical explosions, InSAR, special studies, etc.) or satisfied the Bondár-defined criteria for GT5-25. Posteriori assessments of hypocenter uncertainty and data measurement uncertainty are combined to form datum-specific uncertainty for each tomographic ray. The final data set includes 112,168 distinct events and 2.9 million arrivals. Again, refer to the companion paper above for a more complete description of how the model data set was compiled.

The initial goal for an updated RSTT model is to optimize the results for use with the global IMS network. To that end, we identified phase arrivals specifically from the REB that were associated with the GT events above to use for preferential weighting during the tomography procedures, as well as for later validation steps. The IDC is a component of the Preparatory Commission for the Comprehensive Nuclear-Test- Ban Treaty Organization (CTBTO PrepCom). It receives, collects, processes, analyzes, reports on, and archives the IMS data and distributes IDC products to the State Signatories of the treaty [public access possible via virtual Data Exploitation Centre (vDEC) https://www.ctbto.org/specials/vdec (Vaidya et al. 2009)]. Of the 112,168 tomographic data set events, 7638 of them had at least three regional phases (i.e., Pn, Pg, Sn, or Lg) listed in the REB. To allow for validation events after tomographic inversion as well as data for developing the path-dependent travel-time uncertainty, we randomly selected 30% of the events (2309) to leave out of the tomography data set, resulting in a sufficient number of validation events and regional phases (see Fig. 4 of the tomography companion paper).

3 Path-Dependent Travel-Time Uncertainty

For travel-time tomography models, traditional travel-time path uncertainty or “model” error is typically calculated from a data set that was held out of the tomography process for validation, maintaining non-circularity. These travel-time data are then run through the resulting tomography model and travel-time residuals are calculated and binned into groups by distance between the source and receiver. These residual bins are then used to determine a standard deviation or percentile by distance for estimating the travel-time uncertainty (in seconds) for a given source–receiver path distance (Fig. 1, left).

Fig. 1

(left) Pn travel-time prediction uncertainty using one-dimensional (1D) (with distance) values calculated from the 68^th percentile of travel-time residuals for iasp91 (dotted) and the RSTT validation data (solid). (right) Two-dimensional result using the 1D RSTT values for uncertainty. The 1D, distance-dependent values create a radially-varying pattern from a station, with the same values used for any station

The distance-dependent uncertainty methodology is currently a standard for seismic location (e.g., Kennett and Engdahl 1991). Given the 2.5D model framework for RSTT, we have the means to improve upon the standard. Since the 1D travel-time uncertainty would only vary with distance between the source and receiver, using this uncertainty for a tomography model with spatial variations in 2D or 3D results in the same uncertainty values calculated for any station in any region of the tomography model, no matter the model complexity or data coverage (Fig. 1, right). For a 2D model like RSTT (or 2.5D if the gradient component is considered), we want to have a path-dependent error model with a similar dimension as the model itself. This accounts for spatial variations in empirical data which could lead to bias in the model depending on the specific ray path. Evaluating the variations in residuals among model nodes will allow for a separation of the error into “model”, “bias”, and “random” components, thereby permitting a path-dependent calculation of travel-time uncertainty. The tomographic model is global in extent, but there are many nodes where there is limited or even no ray path coverage. The error model must account for this in some way and convey the lack of data and thus higher uncertainty to the user for any ray paths that traverse these nodes with no data.

Others have calculated the full 3D model covariance matrix (e.g., Ballard et al. 2016a; Hipp et al. 2012) along with 3D tomography, a computationally-efficient spike process to produce a full 3D covariance matrix for large models (Simmons et al. 2019), or more of a brute force data-spike process (one matrix row at a time) (Soldati et al. 2006). These methods work well for estimating path-dependent travel-time uncertainty, but are computationally expensive and are difficult to use in real-time to determine an uncertainty value for a given path. The goal for using RSTT is to provide the general user with a path-dependent, travel-time prediction error on-the-fly. For a 2.5D model at regional distances, we strive to get close to the same prediction uncertainty as when using a full 3D covariance matrix, that is, matching the target percent value for a given coverage ellipse (e.g., 90%).

3.1 Constructing a 2D Random Effects Model for RSTT

A general linear model representation that includes physical bias and error is:

${\text{Observation}}\, = \,Model\, + \,Bias\, + \,RandomError.$

For the development here, this basic model is extended to greater resolution as:

${\text{Observation}}\, = \,Model\, + \,FixedBias\, + \,DynamicBias\, + \,RandomError.$

The Model term, here, is determined external to the equation with geophysical theory (i.e., a geophysical model). However, this geophysical model term will certainly be inadequate in fully describing the travel-time from source to receiver. The term FixedBias is a constant (mean) and assists the Model in predicting travel time or slowness. Intuitively, if the physical Model is good, then FixedBias will be small. The term DynamicBias is defined to account for poorly determined geophysics effects (i.e., model inadequacy) (Anderson et al. 2009) that change from event to event such as local-to-source geology and event depth, and can conceptually be viewed as ModelError. FixedBias is calibrated and ModelError (DynamicBias) is treated as a normal random variable with zero mean. The RandomError term is a combination of arrival pick error as well as other possible errors (e.g., event time, clock drift, localized noise, etc.) and is also normal with zero mean.

If the predictive Model is moved to the left side, the model is: ${\text{Observation-Model}}\, = \,FixedBias\, + \,ModelError\, + \,RandomError,$

and the right side of this equation is a statistical random effects model (REM) (see "Design and Analysis of Experiments", Eighth Edition, Montgomery 2012 for a detailed development of equations). Recall that ModelError and RandomError terms are zero mean so the FixedBias term can be estimated from the REM. As the slowness Model is improved, the various components for FixedBias and ModelError will go down and depending on the improvements, FixedBias could go down more than ModelError. This feature of the REM provides the capability to prioritize research to improve a slowness model (Model).

When we apply the above for the construction of a 2D error model for RSTT, we start with the standard travel-time formulation (Eq. 1)

$$t_{i} = t_{o} + (T_{oi} + \mu_{T} + E_{T} ) + \varepsilon_{i}$$

(1)

where t_i = the phase arrival time, t_o is the event origin time, T_oi is the model-based travel time for event o and arrival i, µ_T is the fixed bias on the travel-time path (T), E_T is the model error for the travel-time path, and ε_i is the random error for the arrival (or also measurement error). Because we are solving for slowness in the tomographic inversion (i.e., inverse velocity), we wish to model the error components for those slowness values instead of travel time. We can use the same general error model above for multiple estimates (observations) of apparent slowness (inverse velocity, [s/km]) local to a node n in the tomography model along a path i,

$$s_{i} = (S_{n} + \mu_{{s_{n} }} + E_{{s_{n} }} ) + \varepsilon_{i}$$

(2)

where s_i is the apparent arrival slowness (Eq. 4 below) of the arrival i over its ray path, S_n is the model slowness assigned to a node n that is touched by the ray path (i.e., has a non-zero path weight), and $\mu_{Sn}$, $E_{Sn}$, and $\varepsilon_{i}$ are analogous to Eq. (1), but for the path slowness values at node n. Moving the model slowness to the left-side of Eq. (2) gives,

$$s_{i} - S_{n} = \mu_{{s_{n} }} + E_{{s_{n} }} + \varepsilon_{i}$$

(3)

where the left side of the equation is now the apparent slowness residual of arrival/path i, touching node n, and the right side of the equation is now the sum of error components (fixed bias, model error, random error).

We attempt to isolate the head-wave portion of the ray path (mantle for Pn/Sn or mid-crust for Pg/Lg) for uncertainty estimates. The apparent observed slowness for the head-wave portion of the ray path is based on the observed travel time of the arrival with estimates of travel-time components removed (i.e., the crust and any upper mantle gradient portion) divided by the head-wave distance (Fig. 2),

$$s_{i} = \frac{{(t_{i} - t_{o} ) - (T_{oi,c} + T_{oi,r} + T_{oi,g} )}}{{d_{i} }}$$

(4)

Fig. 2

Example of ray path through tessellation grid points for a source (o) (star) and a receiver (r) (triangle), displaying tessellation triangles between points (gray) that are “touched” by the head-wave portion of the ray path (i). The total travel time is the sum of the source-side crustal leg (T_oi,c), the receiver-side crustal leg (T_oi,r), the mantle gradient portion (for a Pn or Sn path, not shown), and the head-wave portion (T_oi,h). Only the tessellation nodes touched by the head-wave portion of the ray path are used for calculating the error components

where (t_i-t_o) is the observed travel time, and T_oi,c, T_oi,r, T_oi,g are the calculated travel times for an event o, along ray path i, for the source crustal leg c, the receiver crustal leg r, and the mantle gradient portion g, and d_i is the head-wave distance (km) for path i.

3.2 Compiling Slowness Residuals for Estimating Travel-Time Uncertainty

For every node n touched by the ray path i, the calculated node slowness (S_n) is generalized for the REM to be the calculated apparent slowness of the full ray path (S_i), using the calculated travel time of the head-wave portion (h) of the ray path, divided by the head-wave distance (Fig. 2):

$$S_{i} = \frac{{T_{oi,h} }}{{d_{i} }}.$$

(5)

The final apparent slowness residual is the difference between Eqs. (4) and (5):

$$r_{i} = s_{i} - S_{i} = \frac{{(t_{i} - t_{o} ) - (T_{oi,c} + T_{oi,r} + T_{oi,g} + T_{oi,h} )}}{{d_{i} }}$$

(6)

which is simply the observed travel time minus the calculated travel time, divided by the head-wave distance. For a Pg or Lg ray path, the head-wave distance (d_i) is replaced with the distance along the mid-crustal layer meant to be used to compute travel time for the crustal waveguide, and the gradient portion of the calculated travel time (T_oi,g) is zero. After tomography, the validation data set (the data that were left out of the tomographic inversion) is used to calculate slowness residuals for a given phase that are assigned to all “touched” nodes n of the model along each ray path i:

$$r_{i,n} = s_{i} - S_{i,n}$$

(7)

Using these validation ray paths for an arrival phase, each node n will have a collection/array of apparent slowness residuals assigned to it (Fig. 3) and can have a variable number of ray paths that interact with it; thus, each apparent residual slowness array will have a variable number of values. These collections of apparent slowness residuals will be the basis for calculating the fixed bias, model, and random error components for each node n.

Fig. 3

Example of arrays of apparent slowness residuals (r_i,n) collected for each node n (here numbered from 1 to 7), from any subset of paths i through m. A node can have different numbers of ray paths than other nodes (each m can be different for each node), or can have no ray paths (denoted here as “- “)

3.3 Using the Random Effects Model

We use the above REM on the nodal arrays of apparent slowness residuals to separate the slowness components of error into fixed bias, model error, and random error that will be used to calculate overall travel-time uncertainty for an arrival phase and given ray path. Using the Matlab algorithm “anovan” (Dunn and Clark 1974; Goodnight and Speed 1978; Seber and Lee 2003), we systematically evaluate the arrays of apparent slowness residuals at each node n (Fig. 3). Along with node n, we gather the apparent slowness residuals for nodes that are tessellation neighbors to the central node. Each node is also defined as the center of a “hexagon”, with a node neighbor at the vertex of each side. The group of seven nodes then is also defined as belonging to a hexagon. Both the node itself and the hexagon it is assigned to are defined as random in the anovan algorithm.

When running the REM, it is conceptual to think of the RandomError as variations of the apparent slowness residuals at a node, and the ModelError as variations between the nodes and hexagons. The anovan random effects algorithm returns ModelError and RandomError variance components [s²/km²] for each node grouping, along with a determination of slowness bias [s/km]. These returned values are then assigned to the center node whose position dictated the node neighbors. In order to perform the random effects calculation, a minimum of two nodes in a set of neighbors must have apparent slowness residuals assigned to them.

3.3.1 Building a Grid for the Random Effects Calculation

The assumptions that are built into the RSTT algorithm require that for a given source–receiver distance, the assumed turning-point depth of the ray will depend on the upper mantle velocity and the gradient, and is particularly sensitive to the upper mantle gradient. If the turning-point depth becomes too deep, the RSTT assumption of the constant upper mantle gradient breaks down (Myers et al. 2010). As the gradient effect on the travel times is not part of the random effects calculation above, we chose to subset the data used for determining the error components of the RSTT model into 2° distance bins (from 2° to 16°) after analysis of changes to the error components with larger and smaller distance bin sizes. This permits a focused calculation of the error components for similar-distance rays.

For each grid node in the RSTT model, we select the data that were “touched” by ray paths of a particular distance. Using only those data, we run the REM model and collect the resulting parameter estimates for each node; slowness model error variance τ², random error variance σ², and fixed bias values μ. Each grid node may then have variance and fixed bias values for the full number of distance bins, depending on the distance bin size, or only those distance bins where data are available. After all variances and biases are calculated for each node and distance bin, we extrapolate variance/bias values for those nodes where no data were available. This is performed twice to fill-in variances/biases for a significant portion of possible nodes.

The random effects algorithm used in Matlab anovan can return a ModelError variance component that is negative when minimal data are available for estimation. This issue is rare and, in these cases, the ModelError variance is added to the RandomError variance and the ModelError component is then set to zero. This constrains the model error to be positive, while keeping the total calculated error from the REM (i.e., model + random + bias) the same. Figures 4, 5, 6, and 7 show the slowness model error and bias values of the upper mantle (Pn, Sn) and mid-crustal (Pg, Lg) layers for the 4°, 8°, and 12° distances bins. Model error and fixed bias errors are largest for the shorter distance bins (e.g., 4°) and tend to stabilize by at least the 8° bin, still allowing small fluctuations. The fixed bias values for Sn at the 8° bin (Fig. 6) shows a predominance of negative values suggesting that that the RSTT model does not fit the Sn ray path travel time as well for this set of distance bins. The model would need to have slower upper mantle velocities.

Fig. 4

Maps of model error (left column) and fixed bias (right column) for Pn at 4°, 8°, and 12° distance bins. Error components are calculated at 2° bins from 2° to 16°. Model error is displayed as log₁₀ to better delineate spatial changes. Areas with no data coverage are shown as lighter gray, zero values for model error (see text) are shown as dark gray. Both model error and fixed bias components tend to be highest for shorter distance bins

Fig. 5

Maps of model error (left column) and fixed bias (right column) for Pg at 4°, 8°, and 12° distance bins. Plot features are similar to Fig. 4. Data coverage for Pg is more limited than for other regional phases

Fig. 6

Maps of model error (left column) and fixed bias (right column) for Sn at 4°, 8°, and 12° distance bins. Plot features are similar to Fig. 4. Fixed bias values for the 8° distance bin show large areas of negative values

Fig. 7

Maps of model error (left column) and fixed bias (right column) for Lg at 4°, 8°, and 12° distance bins. Plot features are similar to Fig. 4. Similar to other phases, component values start higher and stabilize to lower values as distance increases

Additionally, we want non-null variances/biases for all nodes at each distance bin, so the final step is to set values for nodes that are null due to a paucity of ray paths. We determine a mean plus two standard deviations at each distance bin for the model variances and for the fixed bias values. These background values are then assigned to the remaining nodes with null values (Fig. 8). Thus, a complete set of surfaces is built for model error, random error, and fixed bias, with the distance bins acting as the z-dimension. Importantly, these approximate solutions to construction of a complete global uncertainty matrix will not be necessary as resolution of our global slowness model is increased with more data.

Fig. 8

Default errors (s/km)² used at nodes for each distance bin where no data are available. Default values are calculated as the mean plus two standard deviations of the error values available for each distance bin

3.4 Calculating a Path-Dependent Travel-Time Uncertainty

Using the complete grid of variances and fixed biases for each distance bin, we can now calculate a full path-dependent travel-time uncertainty estimate for a given ray path (i) of a certain distance between source and receiver (∆). For the ray path (i) (where there are n head-wave nodes that are touched by the path, and m crustal nodes touched by the crustal-leg portions of the path), we gather the nodal weights (in kilometers) from the tomography model:

$$\left[ {\begin{array}{*{20}c} {w_{1} } & {w_{2} } & \cdots & {w_{n} } \\ \end{array} } \right]$$

(8)

$$\left[ {\begin{array}{*{20}c} {\nu_{1} } & {\nu_{2} } & \cdots & {\nu_{m} } \\ \end{array} } \right]$$

(9)

where w are the mantle or mid-crustal head-wave nodal weights [km] and v are the crustal-leg nodal weights [dimensionless].

Using the source–receiver distance (∆), we create a new grid of values for the model variance, the random variance, and the fixed bias as linear-interpolated values from the bounding distance bin grids for each of the variances/bias grids. For example, if the source–receiver distance is 8.65°, and the distance bin is 2°, a new grid for each of the variances/bias is created by interpolating 67.5% of the 8° grids and 32.5% of the 10° grids. These new grids are the only ones necessary for completing the travel-time uncertainty estimate for the current ray path.

The path-based covariance matrix in Eq. (10) is determined using the new interpolated grids above. The general physical-basis structure of the uncertainly matrix in Eq. (10) can be derived from Figs. 2 and 3, and the Appendix in Anderson et al. (2009). Using the notation here, the statistical calculation of mean-squared error can be expanded as estimated FixedBias squared plus the sum of estimated ModelError and RandomError variance components, as given in the diagonal elements of the uncertainty matrix in Eq. (10). For nodes that are touched by the ray path (non-zero node weights), we assume non-zero covariance elements (off-diagonal) between node neighbors only, taking the average of the model error variances for each node. This is an approximate approach to transition between hexagonal regions, and is a topic of continuing research. Given that node spacing is nominally 1°, the constraint of only node-neighbor tomographic smoothing, and that error uncertainty is unlikely to be significant at greater distances, the authors felt that limiting off-diagonal nodal values to be only between neighbors was a valid assumption.

The total travel-time uncertainty is a combination of the uncertainty from the head-wave portion of the ray path (mantle or mid-crust) and the uncertainty from the crustal legs of the ray path. Using a standard double integral calculation (or matrix multiplication), a single estimate of travel-time uncertainty U [s²] for the head-wave portion of a ray path is obtained by multiplying the head-wave nodal weights [km] by the path uncertainty matrix [s²/km²], and again by the head-wave nodal weights [km]. This calculation is formally given in Eq. (10). From this calculation, U is conceptually a matrix mean-squared calculation—a weighted sum of node mean-squared error values and the sum weighted covariance values between these nodes.

Finally, the RandomError variance component [s²], which can be attributed to arrival-time measurement (pick) error, may be included in the final calculation as an estimate of nominal pick uncertainty. If pick uncertainty is added in the location algorithm, then the random error component should be removed from the model-uncertainty calculation so that pick error is not accounted for twice.

$$U_{hw} = \left[ {\begin{array}{*{20}c} {w_{1} } & {w_{2} } & \cdots & {w_{n} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\tau_{1}^{2} + \mu_{1}^{2} + \left( {\sigma_{1,opt}^{2} } \right)} & {\tfrac{{\left( {\tau_{1}^{2} + \tau_{2}^{2} } \right)}}{2}} & \cdots & {\tfrac{{\left( {\tau_{1}^{2} + \tau_{n}^{2} } \right)}}{2}} \\ {\tfrac{{\left( {\tau_{2}^{2} + \tau_{1}^{2} } \right)}}{2}} & {\tau_{2}^{2} + \mu_{2}^{2} + \left( {\sigma_{2,opt}^{2} } \right)} & \cdots & {\tfrac{{\left( {\tau_{2}^{2} + \tau_{n}^{2} } \right)}}{2}} \\ \vdots & \vdots & \ddots & \vdots \\ {\tfrac{{\left( {\tau_{n}^{2} + \tau_{1}^{2} } \right)}}{2}} & {\tfrac{{\left( {\tau_{n}^{2} + \tau_{2}^{2} } \right)}}{2}} & \cdots & {\tau_{n}^{2} + \mu_{n}^{2} + \left( {\sigma_{n,opt}^{2} } \right)} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {w_{1} } \\ {w_{2} } \\ \vdots \\ {w_{n} } \\ \end{array} } \right]$$

(10)

For the crustal-leg portion of the travel-time uncertainty calculation, we currently assume an a priori value using a percent of the calculated source and receiver legs of the crustal travel time, as it is not part of the tomography inversion other than a crustal modifier (Myers et al. 2010). The formulation of the travel-time uncertainty from the crustal legs is similar to Eq. (10) above but there are no off-diagonal components to the variance:

$$U_{crust} = \sum\limits_{j = 1}^{{m_{c} }} {\left( {\nu_{j} p_{j} T_{oi,c} } \right)^{2} } + \sum\limits_{j = 1}^{{m_{r} }} {\left( {\nu_{j} p_{j} T_{oi,r} } \right)^{2} } ,$$

(11)

where v is the dimensionless crustal-leg path weights, p is the percent of the crustal-leg travel time, and T_oi is the crustal leg travel time for an event o and path i for the source crustal leg c and the receiver crustal leg r. Currently, we use a constant value for the percent p, based on a percentage of the crustal-leg travel time (15% for Pn/Pg, 15% Sn/Lg) (e.g., Ballard et al. 2016a; Rodi and Myers 2013). Future values for this a priori percent of crustal travel time could vary by model node, hence the subscript j. For a surface event, the percent of the ray path that traverses the crust will be larger than for an event at depth. The receiver-side crustal leg distance is fairly consistent and only changes with the ray parameter. Thus, events at depth will have a smaller component of travel-time uncertainty from the crust (Fig. 9).

Fig. 9

RMS of the crustal-leg component of travel-time uncertainty for each validation event regional phase relative to the depth traversed in the crust (event-side head-wave (Moho or Mid-Crust) depth—event depth). Negative traversed depth values indicate events below the head-wave interface. Note the reduction of uncertainty as depth reaches the interface depth (i.e., 0 km). For events below the interface, the crustal-leg uncertainty is entirely due (for Pn/Sn) or mostly due (for Pg/Lg) to the receiver-side leg of the ray path

Figure 10 is an example of the resulting components of the travel-time uncertainty calculation, for the Pn phase at IMS station GERES in Germany, and events at zero depth. The resulting total travel-time uncertainty is very different compared to the original IDC source-specific station corrections (IDC-SSSCs) (Firbas et al. 1998) uncertainty values that show a large northeastern area of the surface with low uncertainty. Station GERES is located in an area of the tomography and validation data sets where there are many crossing paths. The total travel-time uncertainty [s] is the square-root of the sum of the squares of the model uncertainty (that has both diagonal and off-diagonal components) and the fixed bias (Eq. 10), along with the crustal uncertainty component for the crustal legs (Eq. 11). For an event at zero depth, the most significant portion of the total uncertainty generally comes from the crustal legs, with the uncertainty value directly related to the crustal thickness at a given event/station location. If an event is much deeper in the crust (e.g., 15 km), the portion of the total uncertainty from the crustal legs is significantly lower (Figs. 9, 11).

Fig. 10

Example of the travel-time uncertainty components (s) for the Pn phase around IMS station GERES and events at zero depth. a Model component, b Fixed Bias component, c Crustal-leg component, d Total Travel-time uncertainty, e Original IDC-SSSC total uncertainty surface (Firbas et al. 1998). Station GERES is located in a part of the tomography/validation data set where there are many ray paths for the various distance bins

Fig. 11

Comparison of Pn crustal leg travel-time uncertainty component at varying event depths for station GERES. (left) 0 km depth, (right) 15 km depth

Owing to the calculation of values at distinct source–receiver distance bins, ring-like patterns may appear depending on the spatial variation of the error components near a station. Figure 12 is an example of the resulting Pn phase travel-time uncertainty components for station BDFB in an area of the model where there are far fewer ray paths in the tomography/validation data set. The component values for BDFB are generally larger than for GERES, with the model component defaulting to the background value for the majority of the 4° distance bin. Figure 13 contains examples of uncertainty components for IMS station GERES and the Pg, Sn, and Lg phases. The model error component of the uncertainty is a more significant portion of the total uncertainty for these phases, when compared to the Pn phase values. This pattern of higher uncertainty values for Pg, Sn, and Lg phases is fairly consistent for all IMS stations.

Fig. 12

Example of the travel-time uncertainty components (s) for the Pn phase around IMS station BDFB and events at zero depth. a Model component, b Fixed Bias component, c Crustal-leg component, d Total Travel-time uncertainty. Station BDFB is located in a part of the tomography/validation data set where there are few ray paths for the various distance bins, especially the 4° bin, where the Model component defaults to the background value for the majority of the distance bin

Fig. 13

Examples of travel-time uncertainty components (s) for the Pg phase (top row), Sn phase (middle row), and Lg phase (bottom row) for IMS station GERES. a Model component, b Fixed Bias component, c Crustal-leg component, d Total Travel-time uncertainty. Events are at zero depth. The Model and Fixed Bias surfaces show higher values overall and are a more significant portion of the total uncertainty when compared to the Pn phase (Fig. 10)

4 Validation of RSTT Model and 2D Path-Dependent Uncertainty

4.1 Validation Arrival Data

We select validation events from a data set developed using the Bayesloc relative relocation algorithm (see “Updates to the Regional Seismic Travel Time (RSTT) Model: 1. Tomography”, this issue) where the events were also listed in the IDC REB and had a minimum of three regional phases (i.e., Pn, Pg, Sn, Lg) within the REB. The event hypocenters were set to be those determined during the Bayesloc relative relocation procedure, while the arrival times were extracted directly from those reported in the REB. All validation events are assumed to be accurate to 5 km or better (i.e., GT5) based on uncertainty estimates after relocation. All regional phases were set to “defining” regardless of whether they were originally defining in the REB.

During the relocation process, the total travel-time uncertainty is a combination of both traditional “model” uncertainty and “pick” uncertainty. ${\text{Total Uncertainty}}^{{2}} \, = \,\left( {"pick"{\text{uncertainty}}} \right)^{{2}} \, + \,\left( {"model"{\text{ uncertainty}}} \right)^{{2}} .$

The Bayesloc-relocated data used for tomography produces values for pick uncertainty based on phase, station, and event factors, which are constrained by posteriori residual statistics (Myers et al. 2007, 2009). Pick and event location uncertainty were used to weight the ray paths during tomographic inversion, thus having a direct impact on the tomography model. In this validation exercise, when using the RSTT model, we compare the 1D uncertainty estimate to the 2D, path-dependent approach described above. However, the “pick” uncertainty used will be that assigned by the IDC in the REB.

Figure 14 shows the root-mean-squared (RMS) values of the total RSTT 2D path-dependent travel-time uncertainty due to components of error (i.e., model error, fixed bias, crustal error) for each event using all available phases. Each arrival associated to an event will have its own values for the error components, depending on the station and event distance relative to the station. Event depth has no direct effect on the model error and fixed bias components of the total travel-time “model” uncertainty. Patterns of the total “model” uncertainty do show higher RMS values in central China, Australia, and parts of Alaska, with lower RMS values especially apparent at subduction zones and along mid-ocean ridges, most likely from the reduced event-leg crustal error values associated with deeper events and thinner crust described above.

Fig. 14

RMS of the total RSTT travel-time “model” uncertainty from components of error (i.e., model error, fixed bias, crustal-leg error) for each of the 2309 validation events and associated regional phases. Values are sorted so higher values are plotted on top. Subduction zones and mid-ocean ridges have lower overall uncertainty, mainly due to reduced event-side crustal-leg error values associated with deeper events and/or thinner crust (Fig. 9)

4.2 Validation Event Relocation

For relocation of the 2309 validation events, we used 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 and ak135), 3D source-specific station corrections (SSSCs), ray-bending through 3D models, and full RSTT models for travel-time prediction. All depths were 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 phases (Pn, Pg, Sn, and Lg), with a limit of 17° (event-station distance) to be consistent with the RSTT assumptions for turning point depths above the 410 km discontinuity. No phases are removed during relocations because of large residuals, and all primary and auxiliary IMS stations are used.

We compare relocation results based on each regional phase (Pn, Pg, Sn, and Lg) and based on all available regional phases to the Bayesloc-defined locations, assuming they are GT5 or better (i.e., epicenters known to within 5 km). For velocity models, we use (1) the iasp91 model (the 1D model used by the IDC), (2) SSSCs produced by the IDC (IDC-SSSCs) (Firbas et al. 1998) (converted to GeoTess format (Ballard et al. 2016b)), (3) the RSTT model developed above (with validation events removed before tomographic inversion) using a distance-dependent model uncertainty (“RSTT-1Du”) (Pn shown in Fig. 1), and (4) RSTT using the 2D, path-dependent uncertainty described above (“RSTT-2Du”). For the IDC-SSSCs, only about one-third of the primary/auxiliary stations have defined SSSCs; the remaining stations will default to the iasp91 model. We evaluate the effects of the path-dependent uncertainty by determining the percent of relocation instances where the 90% coverage ellipse includes the original GT location, relative to azimuthal gap. The GT5 epicenter uncertainty is taken into account when determining if the GT locations are within the relocated error ellipse (i.e., the relocated error ellipse semi-major/minor axes are increased by the square-root of five). In addition, we compare the coverage ellipse size in relation to both the ellipse percentage as well as the 1000 km² ellipse specification in the CTBT. Mislocation results for the iasp91, IDC-SSSCs, and RSTT models/corrections are described in “Updates to the Regional Seismic Travel Time (RSTT) Model: 1. Tomography”, this issue.

4.2.1 Validation Using Pn Phases

Of the 2309 possible REB validation events, 1395 had a minimum of three available arrivals for Pn-only relocations (Fig. 15). Figure 16 shows the coverage ellipse percentage and ellipse area for relocations with the velocity models and uncertainty estimates specified above. The main metric for using the 2D uncertainty estimates for RSTT is to consistently reach the appropriate size for the 90% coverage ellipse—the GT location is within the relocated error ellipse 90% of the time. Because the number of arrivals used for a relocation does not necessarily indicate the azimuthal gap, we chose to display the ellipse percentage relative to azimuthal gap (Fig. 16). The goal is for a consistent 90% value across all azimuthal gaps for any velocity model and uncertainty estimate used, while evaluating the necessary ellipse size to obtain the 90% metric as well as whether the CTBT 1000 km² specification is obtained.

Fig. 15

Validation events with a minimum of three REB arrivals for Pn-only relocations (gray, 1395). Color-filled circles indicate those events whose relocated 90% coverage ellipse contained the GT event epicenter location. Red open circles are those events that did relocate with the model indicated but the relocated coverage ellipse did not contain the GT epicenter location. Models used (bottom-left) (black) iasp91, (bottom-right) (green) IDC-SSSCs, (top-right) (cyan) RSTT with 1D uncertainty, (top-left) (blue) RSTT with 2D path-dependent uncertainty

Fig. 16

Validation results for Pn-only relocations using the velocity models and uncertainties specified. For RSTT, no validation events were used during tomography. (top) Percent of time the relocation 90% coverage ellipse contains the GT location. The 90% coverage ellipse line is shown. Azimuthal gap bins are every 15° and overlap by 1°. The iasp91 model (black) over-estimates the uncertainty, while the IDC-SSSCs (green) and RSTT 1D uncertainty (dashed blue) underestimate the uncertainty. The RSTT 2D uncertainty (solid blue) reaches 91.0%. (middle) The median 90% coverage ellipse area for each model and uncertainty used. The RSTT with 2D uncertainty is closer to the 90% coverage ellipse and has about the same median ellipse sizes as the IDC-SSSCs. The IDC-SSSCs and RSTT with 1D uncertainty hit the CTBT 1000 km² specification (dotted line) at the ~ 110° and ~ 150° azimuthal gap points, respectively, but do not reach the 90% coverage ellipse metric. (bottom) Histogram of counts per azimuthal gap for the iasp91 model to demonstrate numbers of arrivals in each bin

The results for the iasp91 model suggest that the model uncertainty values are too high, overestimating the ellipse percentage. The IDC-SSSCs and RSTT-1Du results suggest those models underestimate the ellipse percentage. The RSTT-2Du reaches an overall value of 91.0%, slightly over the 90% metric, with a relatively “flat” pattern across the azimuthal gap values. The IDC-SSSCs and RSTT-1Du models result in ellipse sizes that reach the CTBT 1000 km² specification (dotted line) at the ~ 110° and ~ 150° azimuthal gap points, respectively, but do not reach the 90% coverage ellipse metric. Neither the iasp91 nor the RSTT-2Du models reach the 1000 km² ellipse size for any azimuthal gap. This is a surprising result and suggests that the use of a single regional phase type for relocation here is not sufficient to overcome the model uncertainties and permit small ellipses sizes that also encompass the GT location.

4.2.2 Validation Using Pg Phases

Of the 2309 possible REB validation events, only 62 had a minimum of three available arrivals for Pg-only relocations (Fig. 17). 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 18 shows the single-event relocation results for coverage ellipse percentage and overall ellipse size. The ellipse percentage results for the Pg-only events show significant fluctuations over azimuthal gap, most likely from the very low number of available events. Results using the iasp91 model show ellipse percentage values at almost the 90% average value, but also the highest ellipse areas. RSTT-2Du is at 83.3% for ellipse percentage with ellipse area at very near the iasp91 level, but with more variation over azimuthal gap. The percentage for RSTT-1Du is much lower at 81.4% and the IDC-SSSCs result is very low at 55.6%. RSTT-2Du results in a significant reduction in mislocation at large azimuthal gaps (> 270°) over RSTT-1Du, presumably from more appropriate model uncertainties being applied to arrival for those events.

Fig. 17

Validation events with a minimum of three REB arrivals for Pg-only relocations (gray, 62). Plot features are similar to Fig. 15

Fig. 18

Validation results for Pg-only relocations using the velocity models and uncertainties specified. Features of plots are similar to Fig. 16 above. The iasp91 model reaches the 90% metric but also with the largest median ellipse area. The RSTT 1D and 2D uncertainty results have a similar ellipse percentage with the 2D uncertainty a few percentage points higher. The IDC-SSSCs significantly underestimate the uncertainty. The large fluctuations in the ellipse percentage values are due to the limited number of Pg-only validation events available. (bottom) The median mislocation values are shown to demonstrate how RSTT with 2D uncertainty reduces mislocation for larger azimuthal gaps

4.2.3 Validation Using Sn Phases

There are 528 REB validation events that had three available arrivals for Sn-only relocations (Fig. 19). Figure 20 shows the single-event relocation results for Sn-only validation events. The Sn-only results for ellipse percentage demonstrate that the 2D uncertainty values used with RSTT meet the 90% coverage ellipse metric but goes slightly over (91.1%), although all the model versions except iasp91 exhibit a decreasing slope with increasing azimuthal gap, along with significant percentage fluctuations. Like for the Pn-only results, iasp91 again overestimates the uncertainties (92.5%) needed to meet the 90% metric, with RSTT-1Du and IDC-SSSCs again underestimating the needed uncertainties (83.8% and 81.3%, respectively). Only RSTT-1Du ever reaches the 1000 km² ellipse area specification. RSTT-1Du does reach the 90% ellipse percent metric, but only at ~ 150° azimuthal gap, with a lower ellipse area value. Presumably, if the RSTT-2Du uncertainty values were lower and relocation results were more precisely reaching the 90% ellipse metric at the 150° gap, the ellipse areas between 1 and 2D uncertainty would be similar.

Fig. 19

Validation events with a minimum of three REB arrivals for Sn-only relocations (gray, 528). Plot features are similar to Fig. 15

Fig. 20

Validation results for Sn-only relocations using the velocity models and uncertainties specified. Features of plots are similar to Fig. 16 above. Like the Pn-only results, the iasp91 model over-estimates the uncertainty needed, while the RSTT 1D uncertainty and IDC-SSSCs significantly underestimate the uncertainty needed to meet the 90% metric. The RSTT 2D uncertainty has an average of 91.1% slightly overshooting the 90% metric. All the results except iasp91 show a downward slope of percentages with increasing azimuthal gap. Only the RSTT 1D uncertainty model ever reaches the 1000 km² ellipse area. The RSTT 1D uncertainty model does reach the 90% ellipse metric, but only at ~ 150° azimuthal gap

4.2.4 Validation Using Lg Phases

There are 434 REB validation events that had three available arrivals for Lg-only relocations (Fig. 21). Figure 22 shows the single-event relocation results for Lg-only validation events. Similar to Pg-only results, ellipse percentage results for Lg-only events show large fluctuations, with significant fall-off of the average percentage values with increasing azimuthal gap. We chose to limit reporting the average values to those with azimuthal gaps less than the 230° point. Given the extreme downward slope for all the model versions, iasp91 overestimates the needed uncertainty (93.4%), with the IDC-SSSCs achieving an average percent very close to the 90% metric (89.8%). Both RSTT versions are lower than the 90% metric, with RSTT-2Du having an average value of 86.7% (below the 230° bounding line). When using 2D uncertainty, RSTT does show an increase in ellipse percent at 265° which is different than the other models that use 1D uncertainty. Like for Pg-only relocations, the use of RSTT-2Du does also show a significant decrease in mislocation at the 265° azimuthal gap point, corresponding to the sudden improvement in ellipse percent, but still below the 90% ellipse metric. Ellipse area for all the models shows a similar relative pattern of increase with increasing azimuthal gap. The RSTT-1Du model does result in ellipse area dipping below the 1000 km² mark at an azimuthal gap of 150° while maintaining the 90% ellipse metric.

Fig. 21

Validation events with a minimum of three REB arrivals for Lg-only relocations (gray, 434). Plot features are similar to Fig. 15

Fig. 22

Validation results for Lg-only relocations using the velocity models and uncertainties specified. Features of plots are similar to Fig. 16 above. The 230° ellipse percentage bounding line was chosen due to the sudden drop in values for all model versions. Like the Pn-only results, the iasp91 model over-estimates the uncertainty needed. In this Lg-only case, the IDC-SSSCs are only slightly below the 90% metric, given the 230° bounding line shown. The RSTT results show the 2D uncertainty version reaching 86.7% with the 1D uncertainty version reaching 86.0%. The RSTT 2D uncertainty model does show an increase in ellipse percentage beyond 265°, unlike the other models. (bottom) Mislocation relative to azimuthal gap shows a significant decrease in mislocation at the same 265° gap point demonstrating a clear improvement using the 2D uncertainty values for RSTT, although not reaching the 90% ellipse metric

5 Validation Using All Regional Phases

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 “typical” regional event location, we allow any combination of regional phases. Figure 23 shows which of the relocated validation events have their 90% coverage ellipse also contain the GT epicenter. 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 proportion of Sn and Lg phases can get as high as ~ 50% (see Fig. 22 from the tomography companion paper for a map showing phase percentages when using all regional phases).

Fig. 23

Validation events with a minimum of three REB arrivals using any/all regional phases (gray). Plot features are similar to Fig. 15

Figure 24 shows the relocation metric results when using any available REB regional phases. For the ellipse percentage test, no model version attained the average 90% metric over the span of azimuthal gap. The iasp91 model did achieve the highest average value of 87.0%. The RSTT-2Du model was able to reach an average of 82.5%. The IDC-SSSCs and RSTT-1Du only reached 74.3% and 70.1%, respectively. There is a slight depression in the percentage values from azimuthal gaps of ~ 150° to ~ 250° in the iasp91 and RSTT models, with the low and high azimuthal gap sections reaching much closer to the 90% metric. Using all regional phases during relocation does reduce ellipse area below the CTBT 1000 km² for each model tested at varying azimuthal gaps, but only iasp91 and RSTT-2Du models approach the 90% ellipse percent metric while crossing the 1000 km² line at azimuthal gaps of 95° and 120°, respectively.

Fig. 24

Validation results for All Regional Phase relocations using the velocity models and uncertainties specified. Features of plots are similar to Fig. 16 above. The ellipse percentage values for all model versions in this case fail to reach the 90% metric, with the iasp91 achieving an average value of 87.0%. The RSTT with 2D uncertainty is able to reach an average of 82.5%. The IDC-SSSCs and RSTT with 1D uncertainty only reach 74.3% and 70.1%, respectively. There is a slight depression in the percentage values from ~ 150° to ~ 250° in the iasp91 and RSTT models, with the low and high azimuthal gap sections reaching much closer to the 90% metric. Using all regional phases during relocation does reduce ellipse area below the 1000 km² for each model tested at varying azimuthal gaps, but only iasp91 and RSTT 2D uncertainty approach the 90% ellipse metric

Given that the IDC-SSSCs exist only for stations in the northern hemisphere (see Fig. 16 from the tomography companion paper), relocation results for error ellipse percentage should be different for events between the northern and southern hemispheres. Using all regional phases, the coverage ellipse percent changes significantly between hemispheres for iasp91 (88.4% northern, 70.2% southern), IDC-SSSCs (74.6% northern, 70.2% southern), and RSTT-1Du (71.5% northern, 53.3% southern), with the IDC-SSSC results for the southern hemisphere being identical to iasp91 due to lack of SSSCs. For RSTT-2Du, the coverage ellipse percent remains consistent at 82.4% for the northern hemisphere and 83.1% for the southern hemisphere, clearly demonstrating the importance of including 2D path-dependent travel-time uncertainty during relocation.

6 Discussion and Conclusions

Using a random effects modeling algorithm, we are able to calculate model error, fixed bias, and random error components at each model node for all phases, or set default background values where there are no data. The validation data are distance-binned at 2° to allow better error modeling using consistent ray path distances. There are assumptions involved in the formulation and calculation of the path-dependent travel-time uncertainties, most of which are made to allow for real-time calculation of the travel-time uncertainty.

We developed the error model with apparent slowness values, where the formulation of the apparent slowness residuals would result in zero error given a “perfect” model (i.e., synthetic tests—zero travel-time residual). The combination of apparent slowness values allows for the specified data grouping and permits the separation into components of error via the REM. The calculation of the path uncertainty does assume that only tessellation node neighbors have off-diagonal components of the path uncertainty matrix. Extending the off-diagonal influence of node neighbors could be an item of further research and testing.

As more data become available for a particular phase, the calculation of the REMs should result in lower absolute fixed bias values as well as reduced model variance. The model variances and fixed bias errors and the resulting component uncertainty estimates are lowest for the Pn phase which has the most data for tomography and calculation of the error model. Model error and fixed bias errors are largest for the shorter distance bins and tend to exhibit less variability by the 8° bin. While not shown, the patterns exhibited for calculated random error also demonstrate similar patterns comparable to those observed for the model variances as well as the fixed bias. The use of the REM provides information on the resulting tomographic model, about which aspects are contributing the most to travel-time prediction uncertainty, allowing for future algorithm adjustment to help improve the model and resulting predictions.

The formulation of the error model also assumes that removing the calculated crustal-leg travel time values from the observed travel time isolates the head-wave portions of either the mantle or mid-crustal path. Because we can only estimate the observed crustal-leg travel times using the calculated times, there will be some crustal-leg uncertainty remaining which would become part of the REM for the head-wave. Without knowledge of the specific portion of the observed travel time resulting from the crustal-leg paths, that portion of the uncertainty cannot be isolated on its own.

Calibration of the error model is necessary given that estimates of arrival “pick” uncertainty are combined with the path-dependent “model” uncertainties to give the total travel-time uncertainty for event location. When determining event location with a 90% coverage ellipse (used by the IDC), we use the 90% ellipse percentage as a metric and did modify the crustal uncertainty a priori values to allow better matching of the 90% metric using the Pn- and Sn-only validation tests. The Pg and Lg model and fixed bias values then had to be multiplied by a factor of five to allow those relocation error ellipse sizes/orientations to reach 90% for those phase-only validation tests. Any phase-specific results not quite matching the 90% ellipse metric appear to combine when evaluating the results using all phases. Users of the model with these REB-centric uncertainties will possibly need to further calibrate the travel-time uncertainty values if the observed “pick” uncertainty values used vary significantly from those used by the IDC in the REB.

The uncertainty attributed to the crustal legs significantly contributes to the overall travel-time uncertainty estimates, especially for the Pn and Sn phases. The overall uncertainties from both Pg and Lg have larger contributions from the model error variances and the fixed biases. The uncertainty from the crustal legs will directly depend on event depth determination as that also dictates the length of the crustal-leg portions of ray paths. Barring the availability of depth phases or other depth controls, seismic locations typically have a trade-off between depth and origin time, making depth one of the least known event parameters. If an event has no depth control and has its depth fixed either too deep or too shallow, it could significantly affect the overall estimate of crustal-leg uncertainty. If an event is a presumed explosion, it would generally be fixed at zero or near-zero depth and would have a better chance of a proper estimate of crustal-leg uncertainty. Earthquakes will have more ambiguity in depth and thus, more ambiguity in the correct crustal-leg travel-time uncertainty values.

Validation of the 2D path-dependent uncertainty method involves testing whether the resulting 90% coverage ellipses after relocation contain the GT epicenter location 90% of the time, as well as the resulting overall error ellipse size. In order to evaluate the changes in error ellipse area, and fairly compare the resulting values, the coverage ellipse percentage metric either has to be met or all models/uncertainty combinations need to have similar coverage ellipse percent values. The RSTT-2Du model generally results in larger error ellipse areas, but these ellipse sizes are required in order to include the GT location and consistently approach or satisfy the 90% coverage ellipse metric. The iasp91 model is consistently above the 90% line for coverage ellipse percent and suggests overestimating of the travel-time uncertainty, while the IDC-SSSCs and RSTT-1Du underestimate the uncertainty.

RSTT-2Du can result in achieving all relocation metrics tested here as well as those in the tomography companion paper: (1) GT location within the coverage ellipse, (2) the smallest error ellipse size, and (3) the lowest mislocation. Figure 25 shows examples of how RSTT-2D effects the error ellipse when RSTT-2Du has the smallest mislocation (two locations—one on land, one at a mid-ocean ridge), another model has the smallest mislocation but RSTT-2Du has the more appropriate error ellipse, and where RSTT-2Du reduces the ellipse size where it should have been much larger to include the GT epicenter.

Fig. 25

Relocation examples demonstrating the effect of using RSTT with 2D path-dependent uncertainty (RSTT-2Du). (white star) Ground-truth (GT) location. (blue with solid error ellipse line) RSTT-2Du, (cyan with blue dashed error ellipse line) RSTT with 1D uncertainty, (green) IDC-SSSCs, (black) iasp91. (top left and right) RSTT-2Du results in the most accurate relocation epicenter as well as reduced coverage ellipse size. (top-right) event along a mid-ocean ridge with IDC-SSSCs producing the same location and error ellipse results as iasp91 due to lack of SSSCs (gray line shows location of ridge). (bottom-left) RSTT-2Du does not result in the most accurate relocation for this event but the error ellipse size and orientation are more appropriate. (bottom-right) None of the models/uncertainties tested result in the GT location being located within the produced error ellipse, with RSTT-2Du resulting in an error ellipse that is far smaller than for the other models that all use 1D distance-dependent uncertainty

The maps showing which relocations did and didn’t have the GT location fall within the error ellipse (Figs. 15, 17, 19, 21, 23) did not indicate a failing of this metric for any of the models in particular regions, displaying more overall patterns consistent with epicenter clusters. Relocations using all regional phases result in the lowest error ellipse sizes compared to phase-specific relocations, with the RSTT-2Du results crossing the CTBT 1000 km² metric line at 150° azimuthal gap while also nearly meeting the 90% coverage ellipse metric. The iasp91 model also reaches the CTBT metric but at a smaller azimuthal gap and requires larger error ellipses. The Pn-only relocation results did not show any model/uncertainty combination crossing the 1000 km² ellipse size metric line. This is surprising, in that Pn is the most predominant regional phase in the validation data set. The RSTT model does reduce travel-time residuals and results in reduced mislocations over other tested models (see tomography companion paper), but mislocation might be further reduced with smaller required error ellipses (perhaps reaching the 1000 km² metric) if the model resolution could be further refined in areas with sufficient data coverage.

During the Pg- and Lg-only relocation tests, the authors noticed that the effect of using RSTT-2Du resulted in a significant reduction in median mislocation values (compared to RSTT-1Du) for events which had an azimuthal gap greater than about 265°. The Pg-only relocations have far more limited numbers of events, but the results are consistent with those for Lg-only relocations. For Pn and Sn relocations (not shown here), the use of RSTT-2Du resulted in almost no change in median mislocation. This suggests that the 2D uncertainties are providing better estimates for model uncertainty, particularly for Pg and Lg, and, thus, the proper weighting of the arrivals during relocation. Events with larger azimuthal gaps tend to be those with fewer phases, so more appropriate uncertainties will significantly affect the resulting epicenter locations for those events.

The RSTT model combined with the 2D path-dependent uncertainty formulation permits the determination of more accurate travel-time uncertainty values, specific to a region of the Earth. The error ellipse sizes are more appropriate and are much closer to a 90% ellipse, thus providing more confidence in the event location results as well as the region-specific error ellipse sizes and orientations.