Subslab flow beneath subduction zones revealed by multiple-layer shear wave splitting inversion

Shear wave splitting is widely used to probe seismic anisotropy, but its depth resolution is limited. Building on the formulation of Silver and Savage (1994), we apply a Bayesian inversion that explicitly accounts for uncertainty in shear-wave polarization orientation (θ) to resolve multilayer anisotropy from splitting parameters. We apply this approach to source-side S and SKS data from the Cocos subduction zone, and to SKS data from station NNA above the South America subduction zone and station SNZO at the southern Hikurangi margin. Subduction in all three regions is shallow to flat, minimizing dip-angle effects on SKS. The inversion yields tightly constrained fast directions for both upper and lower anisotropic layers. Three-layer inversions show progressive rotation with depth consistent with two-layer solutions, but are not required by the data. In the Cocos system, the upper-layer fast direction is unambiguously aligned with Cocos plate motion in the NNR reference frame, consistent with subduction-entrained flow. In contrast, beneath NNA and SNZO, the upper and lower layers exhibit trench-subparallel and trench-normal anisotropy, respectively—opposite in layering sense to the poloidal–toroidal flow structure predicted by dynamic models. If both layers reside in the subslab mantle, the trench-parallel upper layer flow would decouple the deeper mantle from the slab, raising questions about how the apparent subduction-driven flow is maintained at depth. Alternatively, the upper layer may reflect deformation within or above the slab. Possible sources of trench-parallel anisotropy include frozen-in oceanic lithospheric fabrics, trench-parallel mantle-wedge flow, or inherited fabrics related to nearby continental shear zones. These results highlight the complexity of subslab dynamics and demonstrate the value of probabilistic multilayer inversion for interpreting shear-wave splitting.