Abstract Seismic anisotropy provides crucial information on the stress state and geodynamic processes inside the Earth. We develop a novel adjoint-state traveltime tomography method using P-wave traveltime data to simultaneously determine velocity heterogeneity and azimuthal anisotropy of the subsurface. First, an anisotropic eikonal equation is derived to model first-arrival traveltimes in azimuthally anisotropic media. Traveltime tomography is then formulated as an optimization problem constrained by the anisotropic eikonal equation, which is subsequently solved by the adjoint-state method. Ray tracing is not required. Its high accuracy is achieved by solving the anisotropic eikonal equation and the associated adjoint equation with efficient numerical solvers. In addition, an eikonal equation-based earthquake location method for azimuthally anisotropic media is developed to solve the coupled hypocenter-velocity problem. The tomography and earthquake location methods are applied to central California near Parkfield to test their performance in practice. A total of 1,068,850 first P-wave traveltimes clearly maps the velocity heterogeneity and azimuthal anisotropy in the upper and middle crust. The average P-wave velocity model shows a striking velocity contrast across the San Andreas Fault (SAF). In the upper crust, we find structural anisotropy in the SAF zone and stress-induced anisotropy off the SAF zone. In the middle crust, the fast P-wave velocity directions are generally fault-parallel due to the decreased effect of the maximum horizontal compressive stress. In all, the real-data application suggests that the new adjoint-state traveltime tomography method can be reliably used to investigate anisotropic seismic structures.

We discuss the performance of the multiple-grid model parametrization in seismic tomographic inversion. Rather than mapping the velocity perturbation $ Δc(x) = c(x) − c_0(x) $ on only one regular/collocated grid as many previous studies did, we obtain individual Δc(x) models on multiple grids and generate several updated velocity models during one iteration. The average of all the updated velocity models is considered to be the input model of the next iteration. Different grids should partially/fully shift from each other and/or have different grid spacings to form the multiple-grid model parametrization. The efficacy of the multiple-grid model parametrization is demonstrated through the practical example of imaging the P-wave velocity structure along the San Jacinto fault, which is one of the most seismically active faults in California. A series of synthetic recovery examples shows that the multiple-grid model parametrization generally has a better or comparable performance in capturing the heterogeneous subsurface structures than a collocated grid. The root mean square values of the traveltime residuals in the final tomographic models obtained with the multiple-grid model parametrization are smaller than those with collocated grids. Tomographic results reveal strong heterogeneities in the crust along the San Jacinto fault. Significant velocity contrasts are observable across the fault at shallow depths. A low-velocity anomaly dominates the trifurcation area of the San Jacinto fault from the middle crust to the lower crust. Relatively large earthquakes occurred at the boundaries of low-velocity structures but with high-velocity anomalies nearby. All the results suggest that the multiple-grid model parametrization can be a reliable approach in future seismic tomography studies.

In this paper, we propose a wave-equation-based travel-time seismic tomography method with a detailed description of its step-by-step process. First, a linear relationship between the travel-time residual $ Δt = T^{obs}–T^{syn} $ and the relative velocity perturbation $ δc(x)/c(x) $ connected by a finite-frequency travel-time sensitivity kernel $ K(x) $ is theoretically derived using the adjoint method. To accurately calculate the travel-time residual $ Δt $, two automatic arrival-time picking techniques including the envelop energy ratio method and the combined ray and cross-correlation method are then developed to compute the arrival times Tsyn for synthetic seismograms. The arrival times $ T^{obs} $ of observed seismograms are usually determined by manual hand picking in real applications. Travel-time sensitivity kernel $ K(x) $ is constructed by convolving a~forward wavefield $ u(t,x) $ with an adjoint wavefield $ q(t,x) $. The calculations of synthetic seismograms and sensitivity kernels rely on forward modeling. To make it computationally feasible for tomographic problems involving a large number of seismic records, the forward problem is solved in the two-dimensional (2-D) vertical plane passing through the source and the receiver by a high-order central difference method. The final model is parameterized on 3-D regular grid (inversion) nodes with variable spacings, while model values on each 2-D forward modeling node are linearly interpolated by the values at its eight surrounding 3-D inversion grid nodes. Finally, the tomographic inverse problem is formulated as a regularized optimization problem, which can be iteratively solved by either the LSQR solver or a~nonlinear conjugate-gradient method. To provide some insights into future 3-D tomographic inversions, Fréchet kernels for different seismic phases are also demonstrated in this study.

We present a three‐dimensional (3‐D) hybrid method that interfaces the spectral‐element method (SEM) with the frequency‐wave number (FK) technique to model the propagation of teleseismic plane waves beneath seismic arrays. The accuracy of the resulting 3‐D SEM‐FK hybrid method is benchmarked against semianalytical FK solutions for 1‐D models. The accuracy of 2.5‐D modeling based on 2‐D SEM‐FK hybrid method is also investigated through comparisons to this 3‐D hybrid method. Synthetic examples for structural models of the Alaska subduction zone and the central Tibet crust show that this method is capable of accurately capturing interactions between incident plane waves and local heterogeneities. This hybrid method presents an essential tool for the receiver function and scattering imaging community to verify and further improve their techniques. These numerical examples also show the promising future of the 3‐D SEM‐FK hybrid method in high‐resolution regional seismic imaging based on waveform inversions of converted/scattered waves recorded by seismic array.

We demonstrate the feasibility of high-resolution seismic array imaging based on teleseismic recordings using full numerical wave simulations. We develop a hybrid method that interfaces a frequency–wavenumber (FK) calculation, which provides analytical solutions to 1-D layered background models with a spectral-element (SEM) numerical solver to calculate synthetic responses of local media to plane-wave incidence.This hybrid method accurately deals with local heterogeneities and discontinuity undulations, and represents an efficient tool for the forward modelling of teleseismic coda (including converted and scattered) waves. We benchmark the accuracy of the SEM-FK hybrid method against FK solutions for 1-D media. We then compute sensitivity kernels for teleseismic coda waves by interacting the forward teleseismic waves with an adjoint wavefield, produced by injecting coda waves as adjoint sources, based on adjoint techniques. These sensitivity kernels provide the basis for mapping variations in subsurface discontinuities, density and velocity structures through non-linear conjugate-gradient methods. We illustrate various synthetic imaging experiments, including discontinuity characterization, volumetric structural inversion for the crust or subduction zones. These tests show that using pre-conditioners based upon the scaled product of sensitivity kernels for different phases, combining finite-frequency traveltime and waveform inversion, and/or adopting hierarchical inversions from long- to short-period waveforms could reduce the non-linearity of the seismic inverse problem and speed up its convergence. The encouraging results of these synthetic examples suggest that inversion of teleseismic coda phases based on the SEM-FK hybrid method and adjoint techniques is a promising tool for structural imaging beneath dense seismic arrays.