Conventional Pn-wave imaging techniques, typically based on the ray theory, rely on Pn travel times to resolve the 2-D velocity and azimuthal anisotropy structures of the uppermost mantle. However, these approaches neglect the finite-frequency effects of wave propagation, thereby being limited in resolution and accuracy in resolving anisotropy structures. With the advancement of high-performance computing, full-waveform inversion has been increasingly applied to short-period body waves such as Pn. This approach incorporates full-wave finite-frequency sensitivity kernels and anisotropic parameterizations consistent with the elastic wave equation, allowing for high-resolution imaging of upper mantle azimuthal anisotropy. In this study, we perform a series of synthetic experiments using various testing models to qualitatively compare the performance of full-wave Pn inversion with that of conventional Pn travel-time tomography. Our results demonstrate that the waveform-based method improves the resolution and robustness in recovering both velocity and anisotropic parameters, particularly in identifying sharp gradients and deeper features. This study establishes the feasibility and advantages of full-wave Pn inversion for high-resolution azimuthal anisotropy imaging in tectonically complex upper mantle settings.
azimuthal anisotropyfull-waveform inversionPn wavetravel-time tomographyupper mantleThe author(s) declared that financial support was received for this work and/or its publication. This work was supported by the Foundation of National Key Laboratory of Uranium Resources Exploration-Mining and Nuclear Remote Sensing (East China University of Technology) under grants 2025QZ-YZZ-06 and 2024QZ-TD-12, the National Natural Science Foundation of China under grants 42174063, 92155307, 41904046, 41976046, and the Guangdong Provincial Key Laboratory of Geophysical High-resolution Imaging Technology under grant 2022B1212010002.section-at-acceptanceSolid Earth GeophysicsIntroduction
Seismic azimuthal anisotropy in the upper mantle is typically attributed to the lattice-preferred orientation of olivine crystals developed under shear deformation driven by tectonic plate motions (Savage, 1999). The fast directions of P wave azimuthal anisotropy generally align with the maximum shear strain axis. Therefore, imaging P-wave azimuthal anisotropy can provide critical constraints on the deformation context and geodynamic evolution.
Pn tomography has served as a typical and powerful tool for probing P-wave velocity and azimuthal anisotropy structures of the mantle lithosphere (e.g., Hearn, 1996; Al-Lazki et al., 2003; Liang et al., 2004; Seward et al., 2009; Buehler and Shearer, 2010; Lei et al., 2014; Basu and Powell, 2019; Du et al., 2019; He et al., 2019; Andriampenomanana et al., 2020; Illa et al., 2021; He and Lü, 2021; Lu et al., 2022; Du et al., 2022; Takeuchi et al., 2023; Zhang et al., 2023). In regions dominated by simple shear deformation, the fast-wave directions are likely aligned with shear-zone orientations, such as collision zones in the Tibetan Plateau (Pei et al., 2007), the San Andreas Fault system (Hearn, 1996; Buehler and Shearer, 2014), and the Queen Charlotte–Fairweather transform fault system (He and Lü, 2021), etc. Under pure shear conditions, the fast-wave directions generally coincide with the axis of maximum extension, such as the Baikal Rift (He et al., 2019) and the Basin and Range (Hearn, 1996). However, exceptions have also been observed. For example, no clear correlation is found between fast-wave direction and shear-strain orientation along the North Anatolian Fault (Al-Lazki et al., 2004); the fast-wave directions are not aligned with the plate boundary in southern California (Buehler and Shearer, 2014).
The conventional Pn tomography yields lateral variations of azimuthal anisotropy under the assumption that Pn waves propagate laterally in the uppermost mantle along the Moho. This assumption effectively restricts the inversion to two dimensions, thereby limiting the capability of Pn tomography to recover three-dimensional, high-resolution velocity and anisotropic structure. To overcome this challenge, in contrast, 3D techniques of Pn tomography have been recently explored. For instance, Bao and Shen (2020) used full-wave-based Pn inversion to construct a 3D P-wave velocity model beneath the eastern Tibetan Plateau. While their study did not incorporate full-wave inversion of anisotropy, it demonstrates the potential of waveform-based imaging to capture complex anisotropic structures beyond the capability of traditional ray-based techniques. However, a direct comparison of the performance between full-wave Pn inversion and conventional Pn travel-time tomography, particularly in resolving azimuthal anisotropy, is still lacking. Consequently, the relative advantages, limitations, and implications of the two methods remain unclear.
Motivated by this, we conduct synthetic experiments to qualitatively compare the full-wave Pn inversion and conventional Pn travel-time tomography in imaging azimuthal anisotropy of the upper mantle. We show the respective abilities of each method to recover the magnitude and fast wave direction of azimuthal anisotropy as well as the isotropic P-wave velocity. Our results provide theoretical insights and practical guidance for the new generation of imaging azimuthal anisotropy of the upper mantle.
Experimental setup for method comparison
In this study, using a forward solver for the anisotropic elastic wave equation, we generated synthetic data sets from various input models that incorporate varying degrees of heterogeneity in P-wave velocity and azimuthal anisotropy. The same synthetic data are used for both waveform inversion and travel-time tomography, and both inversions share an identical initial model, source–receiver geometry, and spatial discretization scheme to maintain fairness and consistency in the comparison.
We then compare the two methods according to the inverted P-wave velocity perturbation, the fast wave direction, and the magnitude of azimuthal anisotropy, the accuracy of the recovered anomaly boundaries and fast wave directions is considered as the primary evaluation of the method performance. For the full-waveform inversion (FWI) specifically, we consider the reduction in data misfit as an additional measure of effectiveness, because the accuracy of the FWI model can be evaluated from the mismatch between forward-modeled and observed waveforms. Detailed descriptions of the two inversion methods, including parameterization, objective functions, and regularization, are provided in Supplementary Appendix A (for Pn FWI) and Supplementary Appendix B (for Pn travel-time tomography).
Synthetic testsInitial model and data
All synthetic inversion tests in this study are conducted using a common two-layer isotropic background model constructed in spherical coordinates. As shown in Figure 1, the model consists of two layers. The upper layer represents the crust, extending from the surface to a depth of 30 km, with density, P-wave velocity, and shear-wave velocity of 2,600 kg/m3, 6.5 km/s, and 3.7 km/s, respectively. The lower layer represents the upper mantle down to 155 km depth, with density, P-wave velocity, and shear-wave velocity of 3,500 kg/m3, 8.5 km/s, and 4.6 km/s, respectively. The simulation domain spans from −1° to 9° in both longitude and latitude.
Three-dimensional view (a) and plan view (b) of the computational domain. Open circles denote the locations of synthetic sources, and triangles denote the locations of synthetic receivers. The dashed line marks the velocity interface at 30 km depth. The filled markers in (b) denote the source and stations used for computing Pn-wave sensitivity kernels in Figure 2.
Panel (a) is a three-dimensional scatter plot displaying circles and triangles on the top surface of a cube plotted against longitude, latitude, and depth in kilometers, with values ranging from zero to one hundred fifty. Panel (b) is a two-dimensional grid showing circles along the grid edge and triangles filling the inner grid, mapped by longitude and latitude in degrees.
Full-wave sensitivity kernels of Pn-wave phase delays with respect to elastic parameter A, under different frequency bands and epicentral distances. (a,b) Kernels with ∼300 km epicentral distance filtered in the 2–10 s and 4–10 s bands, respectively; (c,d) Same as (a,b) but for ∼600 km epicentral distance. All kernels represent depth cross-sections along the great-circle path between the source and receiver. Distribution of the sources and receivers are shown in Figure 1. (e) Depth profiles of sensitivity kernels at the midpoint along each ray path from (a–d).
Four color gradient contour plots labeled (a) to (d) display values, ranging from negative blue to positive red, against depth (up to one hundred fifty kilometers) and distance, with panels (a) and (b) covering three hundred kilometers, and (c) and (d) six hundred kilometers. All panels share a horizontal dashed line at approximately fifty kilometers depth. Panel titles indicate frequency intervals, two seconds to ten seconds and four seconds to ten seconds. Adjacent to each contour plot is a color scale bar. Panel (e) is a line chart showing depth versus K_A value for all four panels, each line type corresponding to one panel, aiding comparison of depth-dependent trends.
All numerical simulations are based on a collocated-grid finite-difference method in the spherical coordinates (Zhang et al., 2012). The horizontal grid spacing is 0.008° (approximately 0.8 km), and the vertical grid spacing increases gradually with depth, from 0.2 km at the surface to 1.4 km at the bottom. The time length of simulation is 150 s with a time step of 0.02 s. A 12-layer perfectly matched layer is applied to all model boundaries except the free surface. The observation system includes 49 receivers and 32 earthquake sources (Figure 1b), with a constant source depth at 20 km (Figure 1c). The source time function is a Gaussian integral with a half-rise time of 0.25 s. This setting is available to simulate Pn waves with frequency up to 1 Hz.
The phase-delay measurements from the two frequency bands (2–10 s and 4–10 s) were applied to the Pn FWI to invert for the isotropic and anisotropic variations of the upper mantle. Figures 2a–d illustrate the sensitivity kernels of Pn-wave phase difference with respect to the elastic parameter A (Dziewonski and Anderson, 1981), for different epicentral distances and frequency bands. As shown in Figure 2e, increasing the minimum filter period from 2 s to 4 s leads to a deeper maximum sensitivity, indicating that low-frequency signals are more sensitive to deeper structures.
Conventional Pn travel-time tomography inverts the Pn absolute arrival times measured at frequencies higher than 1 Hz, which is inconsistent with the cross-correlation phase delays inverted in the Pn FWI. We adopted the cross-correlation approach proposed by Takeuchi et al. (2023), measuring the travel time residuals between reference and perturbed waveforms after applying a 1–10 s bandpass filter. The validation in the Supplementary Appendix C demonstrates that the cross-correlation measurements closely match ray-theoretical estimates, confirming the reliability of this method for analyzing Pn-wave travel time residuals.
Test 1 – simple model
We first investigated the performance of the two inversion methods in the case of a simple azimuthal anisotropy. We introduce a 222 km × 222 km wide anisotropic anomaly between 30 and 80 km depth in the upper mantle (see Figure 3). The perturbation parameters of this anomaly include: δA=8GPa, δBc=4GPa and δBs=4GPa, corresponding to a Vp anomaly of δVp/Vp=1.57%, an anisotropy magnitude of ξ=2.17% , and a fast-wave direction ϕf=−22.5°.
Model setup for Test 1 (simple structure). (a) Three-dimensional view of the model, with green hollow circles indicating the seismic sources, black triangles representing stations, the dashed line at 30 km depth marking the velocity interface, and the blue cube denoting the high‑velocity anomaly; (b) Horizontal slice of the model, where the blue area is the surface projection of the anomalous body, and the black segments within the anomalous body represent fast‑wave orientations, with angles measured clockwise from true north; (c) Vertical cross‑section, with black segments within the anomalous body showing fast‑wave orientations, where angles are taken clockwise from the upward direction.
Three scientific diagrams labeled (a), (b), and (c) show the geometry of a blue rectangular volume inside a 3D box and its projections with axes for longitude, latitude, and depth in kilometers, surrounded by triangle and circle markers. The text to the right gives geophysical parameters for the blue block: delta Vp over Vp equals one point five seven percent, xi equals two point one seven percent, and phi sub f equals negative twenty two point five degrees.
The FWI typically requires multiple iterations to converge. However, given the relatively small perturbation magnitude in this test, the cross-correlation travel time shifts of the Pn waves satisfy the linearized approximation well (see Equation A2). Hence, we employ the result from the first iteration of Pn FWI to assess the effectiveness of the method in this study.
Based on the regularization-parameter selection procedure described in Supplementary Appendix A.3 (L-curve trade-off analysis between the normalized residual norm and the normalized model norm), we performed grid search over (s, λ, ϵ) and selected the optimal parameters at the maximum-curvature (knee) point of the L-curve. Figure 4 illustrates the selection process for the Pn FWI workflow. The optimal parameters are s=4e−11, λ=1.5e−12, and ϵ=0.7. The travel-time inversion used the same parameter-selection procedure (not shown), yielding s = 1200, λ=50, and ϵ=0.7. Note that the absolute magnitudes of s and λ depend on the scaling of the forward operator and the regularization operators. All subsequent tests used the same parameter-selection procedure.
Determination of the optimal smoothing parameter and damping parameters. The grid search is performed for various values of s between 2e-11 and 6e-11, λ between 5e-13 and 2.5e-12, and ϵ between 0.5 and 0.9, respectively. The tradeoff curves are plotted in terms of the normalized model norm (x-axis) and the normalized residual (y-axis), where the model norm is M=∥δm∥2 and the residual is R=∥Gδm−δτ∥2 (see Supplementary Appendix A.3). The normalized quantities are Mn=M/Mmax and Rn=R/Rmax, where Mmax and Rmax are the maximum values across the grid search. (a) The tradeoff curves (solid lines) according to different combinations of λ and ϵ varying in their search ranges, respectively, where each line denotes one combination. Each star is the approximate maximum curvature point of the corresponding L-shaped tradeoff curve. A similar strategy is used to determine (b)λ and (c)ϵ, respectively, yet the number of tradeoff curves decreases from (a–c) when the number of undetermined parameters decreases.
Three-panel figure showing L-curves for model regularization. Panel (a) shows multiple blue curves with red stars for different smoothing values, (b) shows curves for different damping values, and (c) displays clustered points for different damping ratios. All plots depict normalized residual versus normalized model norm with varying parameterizations.
Figure 5 shows the inversion results for Test 1 using both Pn FWI and Pn travel-time tomography. For quantitative comparison, we use the depth of the maximum vertical gradient in the inversion result to estimate the boundary of the anomaly. In the horizontal slices from the Pn FWI (Figures 5a,b), the recovered Vp anomaly aligns well with the true model. The fast wave directions closely match the true directions, demonstrating that the Pn FWI can effectively capture both velocity and anisotropy features. The slight reduction and smearing of resulting anisotropy magnitude are presumably due to regularization and inadequate convergence at the first iteration. Moreover, the vertical slices (Figures 5c,d) show that the upper boundaries of both velocity and anisotropy anomaly are well recovered by the Pn FWI, while the lower boundaries are deviated by less than 5 km.
Resolved P-wave velocity and azimuthal anisotropy anomalies from Pn FWI and Pn travel-time tomography. (a,b) Horizontal slice at 56 km depth of the resulting model of Pn FWI; (c,d) Vertical sections along latitude 4° of the resulting model of Pn FWI; (e,f) Same as (a,b) except for Pn travel-time tomography. (a,c,e) display P-wave velocity perturbations, while (b,d,f) show anisotropy magnitude and fast wave direction. Red line segments indicate the inverted fast directions, scaled by anisotropy magnitude; black segments denote the true directions. Black dashed rectangles outline the true anomaly boundaries. Each subplot includes a left-side parameter anomaly profile (path indicated in title): red for inversion, blue for the true model, and black dashed lines for the locations of maximum positive and negative gradients, used to assess boundary accuracy.
Six-panel scientific figure showing contour plots and line graphs analyzing seismic properties. Panels (a), (c), and (e) display color-mapped velocity perturbations (δVp/Vp) with blue indicating low values and red indicating high values; corresponding line graphs on the left show variation with latitude or depth. Panels (b), (d), and (f) display azimuthal anisotropy with grayscale backgrounds and vector arrows, alongside line graphs for anisotropy magnitude (ξ). Color and grayscale bars provide percentage scales for the data.
In contrast, the travel-time tomography can only provide horizontal results (Figures 5e,f), which is relatively consistent with the true model for Vp anomaly and fast wave directions. However, the extent of the recovered anisotropy magnitude is noticeably stretched in the north-south direction.
To further evaluate the waveform-fitting of the FWI model, Figures 6b,c compare the synthetic Pn waveforms for four representative stations across the initial, the true, and the inverted models in two frequency bands, respectively. Within the Pn window, the FWI-derived waveforms match those of the true model more closely and are significantly improved upon those of the initial model. Statistically, the variance of the Pn phase differences is reduced from 0.09 s2 (initial model) to 0.05 s2 (FWI model), indicating a notable enhancement in data fitting.
Comparison of synthetic Pn waveforms before and after inversion for Test 1. (a) Red circles and triangles indicate the earthquake and station locations used for waveform comparison. (b,c) show synthetic waveforms in the 2–10 s and 4–10 s bands for the initial (black), true (red solid), and inverted (red dashed) models. The two vertical dashed lines indicate the start and end times of the Pn window used for inversion. The Pn window start time is set to the TauP-predicted Pn arrival time in the starting model (t0). The window length is chosen based on the bandpass period to capture the first-arriving Pn wavelet: for 2–10 s we use t0 to t0 +3 s, and for 4–10 s we use t0 to t0 +5 s.
Figure with three panels. Panel (a) shows a grid map with longitude and latitude axes, populated by circles and triangles, some highlighted in red to identify specific locations labeled 1# to 4#. Panels (b) and (c) each contain four line graphs labeled 1# to 4#, showing black and red waveform plots over time with vertical dashed lines marking intervals, where panel (b) is for 2 seconds to 10 seconds and panel (c) is for 4 seconds to 10 seconds.Test 2 – laterally heterogeneous model
We then investigated the performance of the two inversion methods in the case of laterally variated anisotropy. To reflect such structural complexity, we introduced two laterally adjacent anisotropic anomalies beneath the velocity interface (see Figure 7). Both anomalies have horizontal dimensions of approximately 333 km × 144 km and extend from 30 km to 80 km depth. The perturbation parameters of the high-velocity anomaly (blue block) include: δA=8GPa, δBc=4GPa, and δBs=4GPa, corresponding to a δVp/Vp=1.57%, an anisotropy magnitude of ξ=2.17% , and a fast wave direction of ϕf=–22.5°. The low-velocity anomaly (red block) has opposite parameters: δA=−8GPa, δBc=−4GPa, and δBs=−4GPa, corresponding to δVp/Vp=−1.57%, ξ=2.2%, and ϕf=67.5°.
Model setup for Test 2 (laterally heterogeneous structure). (a) Three-dimensional view of the model, with green hollow circles representing seismic sources, black triangles representing stations, the dashed line at 30 km depth indicating the velocity interface, the blue cube representing the high-velocity anomaly, and the red cube representing the low-velocity anomaly; (b) Horizontal slice of the model, where blue and red areas represent the surface projections of the high-velocity and low-velocity anomalous bodies, respectively, and the black line segments within the anomalous body represent the fast-wave polarization direction, with angles defined as positive clockwise from true north; (c) Vertical cross-section, where the black line segments within the anomalous body represent the fast-wave polarization direction, with angles defined as positive clockwise from the upward direction.
Scientific figure with three panels showing a 3D model (a), a top-down view (b), and a cross-sectional view (c) of subsurface structures using colored blocks; blue and red blocks with black hatching represent different physical properties. Key properties and values for blue-black and red-black regions are listed on the right.
Figures 8a,b show that the Pn FWI method well recovers the lateral boundaries of both anomalies, particularly the inner boundary between them, which closely matches the true model. Although the outer boundaries of the anisotropy magnitude ξ are smeared to some extent, the recovered fast wave directions ϕf align well with the true directions. The vertical slices (Figures 8c,d) further demonstrate that the upper and lower boundaries of the anisotropy magnitude anomalies are resolved well, and the spatial distribution of fast wave directions also agrees with the true model, suggesting that the Pn FWI retains high capability in imaging lateral heterogeneities.
Resolved P-wave velocity and azimuthal anisotropy anomalies from Pn FWI and Pn travel-time tomography for Test 2. (a,b) Horizontal slice at 56 km depth of the resulting model of Pn FWI; (c,d) Vertical sections along latitude 4° of the resulting model of Pn FWI; (e,f) Same as (a,b) except for Pn travel-time tomography. (a,c,e) display P-wave velocity perturbations, while (b,d,f) show anisotropy magnitude and fast wave direction. Red line segments indicate the inverted fast directions, scaled by anisotropy magnitude; black segments denote the true directions. Black dashed rectangles outline the true anomaly boundaries. Each subplot includes a left-side parameter anomaly profile (path indicated in title): red for inversion, blue for the true model, and black dashed lines in (c,d) for the locations of maximum positive and negative gradients, used to assess boundary accuracy.
Six-panel scientific figure compares variations in seismic velocity (δVp/Vp) and azimuthal anisotropy (ξ) across different locations and depths, using color scales for magnitude, latitude, longitude, and depth axes, and vector fields for directionality.
The synthetic Pn waveforms generated from the FWI model (Figure 9) closely match those from the true model in both frequency bands and show significant improvement over the initial model. Statistically, the variance of Pn phase differences is reduced from 0.14 s2 (initial model) to 0.04 s2 (FWI model).
Comparison of synthetic Pn waveforms before and after inversion for Test 2. (a) Synthetic waveforms in the 2–10 s bands, and (b) those in the 4–10 s bands, are shown for the initial (black), true (red solid), and inverted (red dashed) models. The two vertical dashed lines indicate the start and end times of the Pn window used for inversion.
Two sets of line graphs labeled (a) and (b) show time series data from four trials each, comparing black and red lines over time. Vertical dashed lines indicate two specific time points on the x-axis in each plot. In (a), labeled 2s-10s, both lines display multiple oscillations after the first dashed line. In (b), labeled 4s-10s, the graphs show smoother single oscillations that peak after the second dashed line. The x-axes in all plots are labeled “Time(s)” and each subplot is numbered from 1# to 4#.
In contrast, the results from Pn-wave travel-time tomography exhibit noticeable deficiency in imaging lateral heterogeneities (Figures 8e,f). For example, the extent of the recovered model in the longitudinal direction is significantly underestimated, failing to capture the full shape of the anomalies. Meanwhile, the anisotropy magnitude distribution and fast wave direction deviate remarkably from the true model.
Test 3 – vertically layered model
We further investigated the performance of the two inversion methods in the case of vertically layered azimuthal anisotropy. To simulate this scenario, we introduce two vertically stacked anisotropic anomalies beneath the velocity interface (see Figure 10). The shallow anomaly is located between depths of 30–60 km, and the deeper anomaly between 70 and 100 km. Both anomalies have identical dimensions of 222 km × 222 km × 30 km. The perturbation parameters of the shallow low-velocity anomaly include: δA=−8GPa, δBc=−4GPa, and δBs=−4GPa, corresponding to a δVp/Vp=−1.57%, an anisotropy magnitude of ξ=∼2.2%, and a fast-wave direction of ϕf=67.5°. The deeper high-velocity anomaly has opposite perturbations: Vp/Vp=1.57% , ξ=2.17%, and ϕf=−22.5°.
Model setup for Test 3 (vertically layered structure). (a) Three-dimensional view of the model, with green hollow circles representing seismic sources, black triangles representing stations, the dashed line at 30 km depth marking the velocity interface, the blue cube denoting the high-velocity anomaly, and the red cube representing the low-velocity anomaly; (b) Horizontal slice of the model, where the red area is the surface projection of the low-velocity anomalous body, and the black line segments within the anomalous body indicate the fast-wave polarization direction, with angles measured clockwise from true north; (c) Vertical cross-section, where the black line segments within the anomalous body represent the fast-wave polarization direction, with angles taken clockwise from the upward direction.
Three scientific diagrams illustrate the positions of red and blue cuboid anomalies in a 3D grid representing longitude, latitude, and depth in kilometers. Diagram (a) presents a 3D view, diagram (b) shows a top-down map, and diagram (c) depicts a vertical cross-section. Symbols include triangles and circles marking measurement points. Accompanying text specifies blue and red anomaly values for velocity perturbation, anisotropy, and orientation in percentages and degrees.
Figures 11a–d show that the FWI method well recovers the velocity anomalies and anisotropy magnitude boundaries in the horizontal slice at 56 km depth, with the fast wave directions closely matching the true model. In the vertical sections, the two vertically juxtaposed anomalies are clearly distinguished, and the pattern of fast wave directions remains consistent with the true structure. Although the lower boundary of the shallow anomaly is slightly biased, the upper boundary of the deeper anomaly is resolved accurately, which is critical for delineating the transition between different depth layers. Overall, the results demonstrate the effectiveness of Pn FWI in imaging anisotropic layering with high resolution. In terms of waveform fitting, synthetic Pn waveforms from the FWI model show good agreement with those from the true model and perform significantly better than those from the initial model (Figure 12). The variance of the Pn phase differences is reduced from 0.11 s2 (initial model) to 0.04 s2 (FWI model), further demonstrating the effectiveness of the inversion.
Resolved P-wave velocity and azimuthal anisotropy anomalies from Pn FWI and Pn travel-time tomography for Test 3. (a,b) Horizontal slice at 56 km depth of the resulting model of Pn FWI; (c,d) Vertical sections along latitude 4° of the resulting model of Pn FWI; (e,f) Same as (a,b) except for Pn travel-time tomography. (a,c,e) display P-wave velocity perturbations, while (b,d,f) show anisotropy magnitude and fast wave direction. Red line segments indicate the inverted fast directions, scaled by anisotropy magnitude; black segments denote the true directions. Black dashed rectangles outline the true anomaly boundaries. Each subplot includes a left-side parameter anomaly profile (path indicated in title): red for inversion, blue for the true model, and black dashed lines in (c,d) for the locations of maximum positive and negative gradients, used to assess boundary accuracy.
Panel (a) shows a color map and profile graph of δVp/Vp at 56 kilometers depth with values ranging from minus 1 percent to plus 1 percent, and panel (b) displays azimuthal anisotropy with vector arrows and a corresponding profile, both plotted against latitude and longitude. Panel (c) presents a vertical section of δVp/Vp versus depth and latitude with a dipole color pattern, and panel (d) shows azimuthal anisotropy vectors over two zones as a function of depth and latitude with accompanying profiles. Panel (e) provides another color map of δVp/Vp with profiles plotted against longitude and latitude, and panel (f) displays a matching azimuthal anisotropy map with direction vectors and a profile, both ranging up to 2 percent with a grayscale bar.
Comparison of synthetic Pn waveforms before and after inversion for Test 3. (a) Synthetic waveforms in the 2–10 s bands, and (b) those in the 4–10 s bands, are shown for the initial (black), true (red solid), and inverted (red dashed) models. The two vertical dashed lines indicate the start and end times of the Pn window used for inversion.
Figure consisting of eight line graphs arranged in two columns labeled (a) and (b), each with four numbered rows. Column (a) shows 2 seconds to 10 seconds data and column (b) shows 4 seconds to 10 seconds data. Each subplot displays black and red waveforms plotted against time in seconds, with two vertical dashed lines marking intervals. Vertical axes have a range from negative one to positive one, and horizontal axes span approximately nine seconds per row.
In contrast, the Pn travel-time tomography (Figures 12e,f) captures the shallow low-velocity anomaly only, but fails to resolve the deeper high-velocity structure. This limitation is consistent with the physical nature of Pn waves, whose high-frequency energy is primarily concentrated near the top of the upper mantle, making them less sensitive to deeper anomalies. These results further highlight the limitations of ray-based travel-time methods in resolving vertical layering structures.
Discussion
Real seismic observations inevitably contain data noise that affects both travel-time picking and cross-correlation–based phase delay measurements. To demonstrate the effect of data uncertainty on the inversion results qualitatively would be critical for understanding the robustness of the method in future real-data applications. Here we conducted additional synthetic tests by perturbing the observations of Test 3, in which the target model is relatively complex, with zero-mean Gaussian noise added directly to the measurements. The standard deviation of noise is set to 10% of the standard deviation of the noise-free measurements, which is compatible with the posteriori data residual level in this study and with the observed noise level reported in Bao and Shen (2020). The perturbation levels of the Pn waveform cross-correlation delay-time measurements and the travel-time picks are the same. The results (Figure 13) show that the recovered models are close to those of the noise-free experiments, respectively, indicating that both workflows are robust under this noise level and experimental settings.
Inversion performed on synthetic data with 10% added noise. (a,b) Horizontal slice at 56 km depth of the resulting model of Pn FWI; (c,d) Vertical sections along latitude 4° of the resulting model of Pn FWI; (e,f) Same as (a,b) except for Pn travel-time tomography. (a,c,e) display P-wave velocity perturbations, while (b,d,f) show anisotropy magnitude and fast wave direction. Red line segments indicate the inverted fast directions, scaled by anisotropy magnitude; black segments denote the true directions. Black dashed rectangles outline the true anomaly boundaries. Each subplot includes a left-side parameter anomaly profile (path indicated in title): red for inversion, blue for the true model, and black dashed lines in (c,d) for the locations of maximum positive and negative gradients, used to assess boundary accuracy.
Six-panel scientific figure displays seismological simulation results. Panels (a), (c), and (e) show color-coded maps and vertical plots of relative seismic velocity perturbations, with red indicating positive and blue indicating negative values. Panels (b), (d), and (f) depict azimuthal anisotropy with grayscale magnitude and overlaid vectors showing orientation changes, accompanied by corresponding vertical plots. Each pair of panels examines data at specific latitude, longitude, or depth ranges, as described by axes labels. Color bars indicate percentage variations for each mapping.
It has been known that waveform inversion could be affected by nonlinearity when the starting model is inaccurate, leading to a missed windowing or cycle skipping in cross-correlation measurements. However, this risk may be limited for the Pn FWI. First, there are no seismic phases before Pn to contaminate the Pn window. Second, the Pn waveform used to measure the cross-correlation phase delay is typically simple (e.g., Figure 6) and short and thus less likely to be locked onto an incorrect cycle. Three, the Pn waveform is typically insensitive to the complex crustal structures (Lu et al., 2025) not addressed by inaccurate starting models. Similar to traditional Pn traveltime tomography, a rigorous treatment is to add station and event terms into the inversion to absorb possible inaccuracies in the starting model. In contrast, being the forward problem much faster to solve, the travel-time tomography would be easily tackled with fully nonlinear techniques (e.g., Del Piccolo et al., 2024).
The source time function (STF) is prescribed throughout our synthetic tests, whereas uncertainties in the STF and event origin time exist in real-data applications. However, such uncertainties would not hinder either the Pn FWI or travel-time tomography. The data misfit of Pn FWI using cross-correlation–based phase-delays within the Pn window is insensitive to absolute waveform amplitudes when the Pn waveform data with large ambient noise are discarded in inversion. In addition, including event terms into inversion (both Pn FWI and travel-time tomography) can absorb uncertainties in the event origin time, thereby preventing a uniform timing bias from being mapped into Earth structure.
Conclusion
In this study, we compare the performance of Pn FWI and Pn travel-time tomography in imaging heterogeneous structures with azimuthal anisotropy in the upper mantle, based on three representative synthetic models. The results demonstrate that Pn FWI offers significantly higher accuracy and resolution in imaging both P-wave velocity and azimuthal anisotropy. For the simple anomaly model (Test 1), the Pn FWI successfully recovers the spatial distribution of velocity and anisotropy parameters, with fast wave directions closely matching those of the true model. In contrast, although Pn travel-time tomography captures the velocity anomaly to some extent, its ability is limited in accurately resolving the boundaries of anisotropic structures. In the laterally heterogeneous model (Test 2), Pn FWI accurately distinguishes the geometry and fast wave directions of multiple adjacent anomalies, whereas travel-time tomography shows noticeable deviations in both boundary location and fast wave direction recovery. For the vertically layered model (Test 3), Pn FWI clearly resolves the approximately stratified anisotropic structures and depth-dependent fast wave directions, whereas travel-time tomography primarily reflects the shallow anomaly, given relatively short epicentral distances, say, Δ≤10∘.
Overall, Pn FWI is better suited for high-resolution imaging of upper-mantle structures, particularly in regions with pronounced lateral or vertical structural complexity. In contrast, conventional Pn travel-time tomography remains useful as a robust first-order tool when data coverage is extensive but the available reference model is uncertain, and when computational efficiency is a primary concern. The results of this study validate the feasibility and advantages of applying Pn FWI for high-resolution imaging of upper mantle azimuthal anisotropy. Future work may attempt to merge the strategies of the Pn travel-time tomography and FWI, as well as extend the methodology to multi-phase joint inversion and real observations for more comprehensive imaging of regional upper mantle dynamics.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
BL: Formal Analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review and editing. XB: Conceptualization, Formal Analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing. Y-CS: Methodology, Writing – review and editing. WZ: Methodology, Writing – review and editing.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2026.1748581/full#supplementary-material
Edited by:Tariq Alkhalifah, King Abdullah University of Science and Technology, Saudi Arabia
Reviewed by:Andrea Morelli, Istituto Nazionale di Geofisica e Vulcanologia, Italy
Xuewei Bao, Zhejiang University, China
ReferencesAl-LazkiA. I.SandvolE.SeberD.BarazangiM.TurkelliN.MohamadR. (2004). Pn tomographic imaging of mantle lid velocity and anisotropy at the junction of the Arabian, Eurasian and African plates: Pn tomography of the Middle East. Geophys. J. Int.158 (3), 1024–1040. 10.1111/j.1365-246X.2004.02355.xAl‐LazkiA. I.SeberD.SandvolE.TurkelliN.MohamadR.BarazangiM. (2003). Tomographic Pn velocity and anisotropy structure beneath the Anatolian plateau (eastern Turkey) and the surrounding regions. Geophys. Res. Lett.30 (24). 10.1029/2003GL017391AndriampenomananaF.NybladeA. A.WysessionM. E.DurrheimR. J.TilmannF.BarruolG. (2020). Seismic velocity and anisotropy of the uppermost mantle beneath Madagascar from Pn tomography. Geophys. J. Int.224 (1), 290–305. 10.1093/gji/ggaa458BaoX.ShenY. (2020). Early‐stage lithospheric foundering beneath the Eastern Tibetan Plateau revealed by full‐wave P n tomography. Geophys. Res. Lett.47 (8), e2019GL086469. 10.1029/2019GL086469BasuU.PowellC. (2019). Pn tomography and Anisotropy study of the central United States. J. Geophys. Res.: Solid Earth124 (7), 7105–7119. 10.1029/2018JB016538BuehlerJ. S.ShearerP. M. (2010). Pn tomography of the Western United States using USArray. J. Geophys. Res.: Solid Earth115 (B9). 10.1029/2009JB006874BuehlerJ. S.ShearerP. M. (2014). Anisotropy and Vp/Vs in the uppermost mantle beneath the Western United States from joint analysis of Pn and Sn phases. J. Geophys. Res.: Solid Earth119 (2), 1200–1219. 10.1002/2013JB010559Del PiccoloG.VanderBeekB. P.FaccendaM.MorelliA.ByrnesJ. S. (2024). Imaging upper‐mantle anisotropy with transdimensional Bayesian Monte Carlo sampling. Bull. Seismol. Soc. Am.114 (3), 1214–1226. 10.1785/0120230233DuG.WuQ.ZhangX.HeJ.ZouL.FengY. (2019). Pn wave velocity and anisotropy underneath the central segment of the North-South Seismic belt in China. J. Asian Earth Sci.184, 103941. 10.1016/j.jseaes.2019.103941DuM.LeiJ.ZhaoD.LuH. (2022). Pn anisotropic tomography of northeast Asia: new insight into subduction dynamics and volcanism. J. Geophys. Res.: Solid Earth127 (1), e2021JB023080. 10.1029/2021JB023080DziewonskiA. M.AndersonD. L. (1981). Preliminary reference Earth model. Phys. Earth Planet. Inter.25 (4), 297–356. 10.1016/0031-9201(81)90046-7HeY.LüY. (2021). Anisotropic Pn tomography of Alaska and adjacent regions. J. Geophys. Res.: Solid Earth126 (11). 10.1029/2021JB022220HeJ.LiY.SandvolE.WuQ.DuG.ZhangR. (2019). Tomographic Pn velocity and anisotropy structure in Mongolia and the adjacent regions. J. Geophys. Res.: Solid Earth124 (4), 3662–3679. 10.1029/2018jb016440HearnT. M. (1996). Anisotropic Pn tomography in the western United States. J. Geophys. Res.: Solid Earth101 (B4), 8403–8414. 10.1029/96JB00114IllaB.ReshmaK. S.KumarP.SrinageshD.HaldarC.KumarS. (2021). Pn tomography and anisotropic study of the Indian shield and the adjacent regions. Tectonophysics813, 228932. 10.1016/j.tecto.2021.228932LeiJ.LiY.XieF.TengJ.ZhangG.SunC. (2014). Pn anisotropic tomography and dynamics under eastern Tibetan plateau. J. Geophys. Res.: Solid Earth119 (3), 2174–2198. 10.1002/2013JB010847LiangC.SongX.HuangJ. (2004). Tomographic inversion of Pn travel times in China. J. Geophys. Res.: Solid Earth109 (B11). 10.1029/2003JB002789LuH.LeiJ.ZhaoD.XuY.SunC.HuX. (2022). Pn anisotropic tomography of Hainan island and surrounding areas: new insights into the Hainan mantle plume. J. Geophys. Res.: Solid Earth127 (6), e2021JB023609. 10.1029/2021JB023609LuB.BaoX.SunY. C.ZhangW. (2025). On the strategy of full-wave Pn inversion for upper-mantle azimuthal anisotropy. Geophys. J. Int.242 (1), ggaf177. 10.1093/gji/ggaf177PeiS.ZhaoJ.SunY.XuZ.WangS.LiuH. (2007). Upper mantle seismic velocities and anisotropy in China determined through Pn and Sn tomography. J. Geophys. Res.: Solid Earth112 (B5). 10.1029/2006JB004409SavageM. K. (1999). Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting?Rev. Geophys.37 (1), 65–106. 10.1029/98RG02075SewardA. M.HendersonC. M.SmithE. G. C. (2009). Models of the upper mantle beneath the central North Island, New Zealand, from speeds and anisotropy of subhorizontal P waves (Pn). J. Geophys. Res.: Solid Earth114 (B1). 10.1029/2008JB005805TakeuchiN.IsseT.KawakatsuH.ShiobaraH.SugiokaH.ItoA. (2023). Olivine fabrics in the Oceanic lithosphere constrained by Pn azimuthal anisotropy. J. Geophys. Res.: Solid Earth128 (6), e2023JB026428. 10.1029/2023JB026428ZhangW.ShenY.ZhaoL. (2012). Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method: 3-D synthetics in spherical coordinates. Geophys. J. Int.188 (3), 1359–1381. 10.1111/j.1365-246X.2011.05331.xZhangZ.IrvingJ. C.SimonsF. J.AlkhalifahT. (2023). Seismic evidence for a 1000 km mantle discontinuity under the Pacific. Nat. Commun.14 (1), 1714. 10.1038/s41467-023-37067-x36973245