Comparative analysis of Pn full-waveform inversion and travel-time tomography in imaging upper mantle azimuthal anisotropy

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.

1 Introduction

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.

2 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).

3 Synthetic tests

3.1 Initial 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/m³, 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/m³, 8.5 km/s, and 4.6 km/s, respectively. The simulation domain spans from −1° to 9° in both longitude and latitude.

Figure 1

Figure 1. 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.

Figure 2

Figure 2. 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).

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.

3.2 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: $\delta A=8GPa$, $\delta B_c=4GPa$ and ${\delta B}_s=4GPa$, corresponding to a ${\mathrm{V}}_{\mathrm{p}}$ anomaly of $\delta {\mathrm{V}}_{\mathrm{p}}/{\mathrm{V}}_{\mathrm{p}}=1.57\%$, an anisotropy magnitude of $\xi =2.17\%$ , and a fast-wave direction ${\varphi }_f=-22.5°$.

Figure 3

Figure 3. 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.

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$, $\lambda$, $ϵ$) 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=4\mathrm{e}-11$, $\lambda =1.5\mathrm{e}-12$, and $ϵ=0.7$. The travel-time inversion used the same parameter-selection procedure (not shown), yielding s = 1200, $\lambda =50$, and $ϵ=0.7$. Note that the absolute magnitudes of $s$ and $\lambda$ depend on the scaling of the forward operator and the regularization operators. All subsequent tests used the same parameter-selection procedure.

Figure 4

Figure 4. 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, $\lambda$ 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=\parallel \delta \mathit{m}{\parallel }_2$ and the residual is $R=\parallel \mathit{G}\delta \mathit{m}-\delta \mathit{\tau }{\parallel }_2$ (see Supplementary Appendix A.3). The normalized quantities are $M_n=M/M_{\max }$ and $R_n=R/R_{\max }$, where $M_{\max }$ and $R_{\max }$ are the maximum values across the grid search. (a) The tradeoff curves (solid lines) according to different combinations of $\lambda$ 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) $\lambda$ and (c) $ϵ$, respectively, yet the number of tradeoff curves decreases from (a–c) when the number of undetermined parameters decreases.

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 ${\mathrm{V}}_{\mathrm{p}}$ 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.

Figure 5

Figure 5. 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.

In contrast, the travel-time tomography can only provide horizontal results (Figures 5e,f), which is relatively consistent with the true model for ${\mathrm{V}}_{\mathrm{p}}$ 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 s² (initial model) to 0.05 s² (FWI model), indicating a notable enhancement in data fitting.

Figure 6

Figure 6. 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 ($t_0$). The window length is chosen based on the bandpass period to capture the first-arriving Pn wavelet: for 2–10 s we use $t_0$ to $t_0$ +3 s, and for 4–10 s we use $t_0$ to $t_0$ +5 s.

3.3 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: $\delta A=8GPa$, $\delta B_c=4GPa$, and $\delta B_s=4GPa$, corresponding to a $\delta {\mathrm{V}}_{\mathrm{p}}/{\mathrm{V}}_{\mathrm{p}}=1.57\%$, an anisotropy magnitude of $\xi =2.17\%$ , and a fast wave direction of ${\varphi }_f=$–22.5°. The low-velocity anomaly (red block) has opposite parameters: $\delta A=-8GPa$, $\delta B_c=-4GPa$, and $\delta B_s=-4GPa$, corresponding to $\delta {\mathrm{V}}_{\mathrm{p}}/{\mathrm{V}}_{\mathrm{p}}=-1.57\%$, $\xi =2.2\%$, and ${\varphi }_f=67.5°$.

Figure 7

Figure 7. 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.

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 $\xi$ are smeared to some extent, the recovered fast wave directions ${\varphi }_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.

Figure 8

Figure 8. 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.

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 s² (initial model) to 0.04 s² (FWI model).

Figure 9

Figure 9. 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.

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.

3.4 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: $\delta A=-8GPa$, $\delta B_c=-4GPa$, and $\delta B_s=-4GPa$, corresponding to a $\delta {\mathrm{V}}_{\mathrm{p}}/{\mathrm{V}}_{\mathrm{p}}=-1.57\%$, an anisotropy magnitude of $\xi =\sim 2.2\%$, and a fast-wave direction of ${\varphi }_f=67.5°$. The deeper high-velocity anomaly has opposite perturbations: ${\mathrm{V}}_{\mathrm{p}}/{\mathrm{V}}_{\mathrm{p}}=1.57\%$ , $\xi =2.17\%$, and ${\varphi }_f=-22.5°.$

Figure 10

Figure 10. 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.

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 s² (initial model) to 0.04 s² (FWI model), further demonstrating the effectiveness of the inversion.

Figure 11

Figure 11. 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.

Figure 12

Figure 12. 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.

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.

4 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.

Figure 13

Figure 13. 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.

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.

5 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, $\Delta \le {10}^{\circ }$.

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.

Funding

The 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.

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

References

Al-Lazki, A. I., Sandvol, E., Seber, D., Barazangi, M., Turkelli, N., and Mohamad, R. (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. doi:10.1111/j.1365-246X.2004.02355.x

CrossRef Full Text | Google Scholar

Al-Lazki, A. I., Seber, D., Sandvol, E., Turkelli, N., Mohamad, R., and Barazangi, M. (2003). Tomographic Pn velocity and anisotropy structure beneath the Anatolian plateau (eastern Turkey) and the surrounding regions. Geophys. Res. Lett. 30 (24). doi:10.1029/2003GL017391

CrossRef Full Text | Google Scholar

Andriampenomanana, F., Nyblade, A. A., Wysession, M. E., Durrheim, R. J., Tilmann, F., Barruol, G., et al. (2020). Seismic velocity and anisotropy of the uppermost mantle beneath Madagascar from Pn tomography. Geophys. J. Int. 224 (1), 290–305. doi:10.1093/gji/ggaa458

CrossRef Full Text | Google Scholar

Bao, X., and Shen, Y. (2020). Early-stage lithospheric foundering beneath the Eastern Tibetan Plateau revealed by full-wave P n tomography. Geophys. Res. Lett. 47 (8), e2019GL086469. doi:10.1029/2019GL086469

CrossRef Full Text | Google Scholar

Basu, U., and Powell, C. (2019). Pn tomography and Anisotropy study of the central United States. J. Geophys. Res.: Solid Earth 124 (7), 7105–7119. doi:10.1029/2018JB016538

CrossRef Full Text | Google Scholar

Buehler, J. S., and Shearer, P. M. (2010). Pn tomography of the Western United States using USArray. J. Geophys. Res.: Solid Earth 115 (B9). doi:10.1029/2009JB006874

CrossRef Full Text | Google Scholar

Buehler, J. S., and Shearer, P. 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 Earth 119 (2), 1200–1219. doi:10.1002/2013JB010559

CrossRef Full Text | Google Scholar

Del Piccolo, G., VanderBeek, B. P., Faccenda, M., Morelli, A., and Byrnes, J. S. (2024). Imaging upper-mantle anisotropy with transdimensional Bayesian Monte Carlo sampling. Bull. Seismol. Soc. Am. 114 (3), 1214–1226. doi:10.1785/0120230233

CrossRef Full Text | Google Scholar

Du, G., Wu, Q., Zhang, X., He, J., Zou, L., Feng, Y., et al. (2019). Pn wave velocity and anisotropy underneath the central segment of the North-South Seismic belt in China. J. Asian Earth Sci. 184, 103941. doi:10.1016/j.jseaes.2019.103941

CrossRef Full Text | Google Scholar

Du, M., Lei, J., Zhao, D., and Lu, H. (2022). Pn anisotropic tomography of northeast Asia: new insight into subduction dynamics and volcanism. J. Geophys. Res.: Solid Earth 127 (1), e2021JB023080. doi:10.1029/2021JB023080

CrossRef Full Text | Google Scholar

Dziewonski, A. M., and Anderson, D. L. (1981). Preliminary reference Earth model. Phys. Earth Planet. Inter. 25 (4), 297–356. doi:10.1016/0031-9201(81)90046-7

CrossRef Full Text | Google Scholar

He, Y., and Lü, Y. (2021). Anisotropic Pn tomography of Alaska and adjacent regions. J. Geophys. Res.: Solid Earth 126 (11). doi:10.1029/2021JB022220

CrossRef Full Text | Google Scholar

He, J., Li, Y., Sandvol, E., Wu, Q., Du, G., Zhang, R., et al. (2019). Tomographic Pn velocity and anisotropy structure in Mongolia and the adjacent regions. J. Geophys. Res.: Solid Earth 124 (4), 3662–3679. doi:10.1029/2018jb016440

CrossRef Full Text | Google Scholar

Hearn, T. M. (1996). Anisotropic Pn tomography in the western United States. J. Geophys. Res.: Solid Earth 101 (B4), 8403–8414. doi:10.1029/96JB00114

CrossRef Full Text | Google Scholar

Illa, B., Reshma, K. S., Kumar, P., Srinagesh, D., Haldar, C., Kumar, S., et al. (2021). Pn tomography and anisotropic study of the Indian shield and the adjacent regions. Tectonophysics 813, 228932. doi:10.1016/j.tecto.2021.228932

CrossRef Full Text | Google Scholar

Lei, J., Li, Y., Xie, F., Teng, J., Zhang, G., Sun, C., et al. (2014). Pn anisotropic tomography and dynamics under eastern Tibetan plateau. J. Geophys. Res.: Solid Earth 119 (3), 2174–2198. doi:10.1002/2013JB010847

CrossRef Full Text | Google Scholar

Liang, C., Song, X., and Huang, J. (2004). Tomographic inversion of Pn travel times in China. J. Geophys. Res.: Solid Earth 109 (B11). doi:10.1029/2003JB002789

CrossRef Full Text | Google Scholar

Lu, H., Lei, J., Zhao, D., Xu, Y., Sun, C., and Hu, X. (2022). Pn anisotropic tomography of Hainan island and surrounding areas: new insights into the Hainan mantle plume. J. Geophys. Res.: Solid Earth 127 (6), e2021JB023609. doi:10.1029/2021JB023609

CrossRef Full Text | Google Scholar

Lu, B., Bao, X., Sun, Y. C., and Zhang, W. (2025). On the strategy of full-wave Pn inversion for upper-mantle azimuthal anisotropy. Geophys. J. Int. 242 (1), ggaf177. doi:10.1093/gji/ggaf177

CrossRef Full Text | Google Scholar

Pei, S., Zhao, J., Sun, Y., Xu, Z., Wang, S., Liu, H., et al. (2007). Upper mantle seismic velocities and anisotropy in China determined through Pn and Sn tomography. J. Geophys. Res.: Solid Earth 112 (B5). doi:10.1029/2006JB004409

CrossRef Full Text | Google Scholar

Savage, M. K. (1999). Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting? Rev. Geophys. 37 (1), 65–106. doi:10.1029/98RG02075

CrossRef Full Text | Google Scholar

Seward, A. M., Henderson, C. M., and Smith, E. 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 Earth 114 (B1). doi:10.1029/2008JB005805

CrossRef Full Text | Google Scholar

Takeuchi, N., Isse, T., Kawakatsu, H., Shiobara, H., Sugioka, H., Ito, A., et al. (2023). Olivine fabrics in the Oceanic lithosphere constrained by Pn azimuthal anisotropy. J. Geophys. Res.: Solid Earth 128 (6), e2023JB026428. doi:10.1029/2023JB026428

CrossRef Full Text | Google Scholar

Zhang, W., Shen, Y., and Zhao, L. (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. doi:10.1111/j.1365-246X.2011.05331.x

CrossRef Full Text | Google Scholar

Zhang, Z., Irving, J. C., Simons, F. J., and Alkhalifah, T. (2023). Seismic evidence for a 1000 km mantle discontinuity under the Pacific. Nat. Commun. 14 (1), 1714. doi:10.1038/s41467-023-37067-x

PubMed Abstract | CrossRef Full Text | Google Scholar