The impact of lithospheric rheology and subduction dynamics on intraslab stress and earthquake characteristics, as well as the cause of the scarcity of continental mantle earthquakes remain debated. We investigate the 2006 Pingtung offshore earthquake doublet (Mw 6.9) in the northern Manila subduction zone, where thinned continental lithosphere is subducting. These events exhibit complex source processes at mantle depths, making their fault geometry poorly constrained. Using potency density tensor inversion on teleseismic waveforms, we resolve the source process without predefined fault geometry. The first event exhibits normal faulting with minor reverse faulting at a deeper depth, while the second shows deeper strike-slip faulting. These diverse focal mechanisms reveal scattered mantle faulting and intraslab stress heterogeneity. Geodynamic modeling indicates mantle earthquakes require both strong subducting mantle lithosphere and substantial bending-induced differential stress. The global scarcity of such events stems from insufficient differential stress in most continental lithospheres.
Earthquakes within the continental mantle lithosphere are rare. The reasons for their scarcity remain a topic of ongoing debate. A central question pertains to the strength of the continental lithospheric mantle. The strength of the continental mantle lithosphere is uncertain due to the complex compositional and thermal heterogeneities associated with the prolonged geological evolution of the continents. Some researchers attribute the rarity of mantle earthquakes to the weak rheology, i.e. ductility, of the continental mantle. Others, however, argue that the mantle lithosphere is strong, but the tectonic stress is commonly insufficient to induce earthquakes. Despite this, the sporadic occurrence of continental mantle earthquakes in some regions suggests that the mantle in those regions is strong enough to accumulate the elastic stress necessary for rupture. Additionally, a strong mantle is proposed to be necessary to support the structural evolution of orogenic and subduction zones. However, because continental mantle earthquakes are rare, resolving this debate remains difficult. The northern Manila subduction zone is notable for relatively frequent mantle lithosphere earthquakes (Fig. 1), and represents a rare case where a thinned Eurasian continental crust and lithosphere is subducting under the Philippian Sea plate (Fig. 1). The continental crust in this region is thinned, ranging from 10 to 30 km, due to its earlier tectonic history as the passive margin of the South China Sea. In some areas, the crust is even thinner, with thicknesses less than 10 km. Although the continental lithosphere in this region has undergone thinning, we assume it retains the compositional characteristics of typical continental lithosphere. Therefore, this region serves as a valuable natural laboratory for investigating the subduction of continental lithosphere and the findings from this study may have broader implications for understanding the behavior and strength of continental lithosphere in general.
Two major seismic events occurred at UTC 12:26 and 12:34 on December 26, 2006. These two events are an earthquake doublet, each with magnitudes of 7.1 and 6.9. This doublet is referred to as the 2006 Pingtung offshore earthquake doublet, named after the coastal area in Taiwan. This doublet is not located on the subduction interface but within the subducting slab and composes intraslab events. Such events -- like many intraslab earthquakes -- are potentially more hazardous than typical interplate or crustal faulting events due to generally higher stress drops, greater radiated energy, and richer high-frequency content. They, despite their potential to produce strong ground shaking, are often underrepresented in regional seismic hazard models. These source characteristics, along with their underlying physical controls, remain poorly understood, underscoring the need for deeper insight into intraslab rupture processes and their tectonic context. Due to poor station coverage, uncertainty in seismic velocity models, differences in analysis methods, and the complex rupture process, different studies report different depth estimates on this doublet (Table S1). However, a consistent pattern is observed across different catalogs, with the first event occurring at a shallower depth and the second event deeper (Table S1). The Global CMT project lists the centroid depths at 19.6 km and 32.8 km for the two events, while the USGS catalog places them at 25.5 km and 32.8 km, respectively. Relocated hypocenters from the Taiwan's Central Weather Administration are 44.1 km and 50.2 km, respectively. Depth phase analysis indicates that the hypocenters and largest asperities of this doublet occurred at depths approximately 38 to 56 km, with the first event at a shallower depth than the second event. A recent study combining regional data put the centroid depths at 28 and 44 km. In addition to their separated depths, this earthquake doublet exhibited distinct focal mechanisms. The first event displayed predominantly normal faulting, while the second, deeper event showed strike-slip faulting (Fig. S1). This suggests considerable stress heterogeneity over a short vertical range. Additionally, the doublet's rupture processes, as evidenced by non-double couple components (Fig. S1) -- ~15% and ~68%, respectively -- should highlight the complexity of this event. Previous studies have similarly indicated complex seismic source processes, but no consensus on a kinematic source model has been reached. Given the sensitivity of intraslab stress regimes on the lithospheric rheology and subduction zone dynamics, constructing a comprehensive seismic source model could offer insights into the stress state and rheological parameters of the continental mantle lithosphere and the associated geodynamic implications.
In this study, we first explore seismic source characteristics by constructing a kinematic seismic source model, emphasizing on depth-dependent variations in source mechanisms. We then establish a viable geodynamic model of the seismic source area, whose stress regimes are consistent with that inferred from the constructed seismic source model. The results offer insights into regional lithospheric strength, providing explanations for the occurrence of mantle-depth earthquakes in this region, and aim to shed light on the broader question of why mantle earthquakes are rare in most continental lithospheres.
Due to the complex source processes and no fault geometry derived from surface expressions, developing a satisfactory source model for this doublet has been challenging. Previous studies of source characteristics have not reached a consensus, for the strikes and dips of the faults, and finite fault slip distributions for the doublet. To solve such complex source problem, most of the source studies would first have to construct a detailed fault model, inferred from near source surface displacement. However, this doublet occurred offshore, and there is no observation on surface ruptures. In addition, this area has low background seismicity. Thus, we have no prior knowledge of the fault geometry. Since the Green's function is sensitive to the details of the assumed fault plane, the solution obtained would be highly dependent on the assumed fault plane. An imprecise fault model can potentially compromise the accuracy of the seismic source solution and bias our interpretation. Another challenge is that although there are near-field strong motion stations, the events are outside the network area and have poor azimuthal coverage.
To address this issue, we apply an advanced finite-fault inversion method-Potency Density Tensor Inversion (PDTI)-which estimates the potency density tensor projected onto the assumed model plane, as described in Shimizu et al.. PDTI is a high-degree-of-freedom source inversion technique that uses teleseismic body waves, which offer broad azimuthal coverage and are relatively insensitive to the precise slip location. It represents slip direction on the model plane as a focal mechanism expressed by a linear combination of five basis double-couple components. Importantly, PDTI allows estimation of the slip direction independently of the fixed orientation of the model plane. This model domain can accommodate multiple faulting episodes of an event, without implying a single fault plane extending in a continuous rupture. The method estimates the spatiotemporal distribution of the potency-rate density tensor, including information on slip vectors and fault geometry or any focal mechanism required by the data, within the defined model domain. The strike of each model plane was determined based on focal mechanisms derived from the W-phase solutions in the USGS earthquake catalog. For the first event (12:26 UTC), the model plane extends along a 171° strike, spanning a depth range of 10 to 60 km and a horizontal distance of 60 km, covering the potential depth range for the source processes. For the second event (12:34 UTC), the model plane strikes along 151° and spans a depth range from 20 to 70 km, extending 60 km along strike. Both strike directions are nearly parallel to the Manila Trench, likely reflecting structures associated with slab bending and unbending. For more details, see the Methods section #1.
We then set up 2D thermo-mechanical models focused on the seismic source region, oriented in the east-west direction, using elasto-visco-plastic rheology. Our 2D numerical code, geoflac, employs the Fast Lagrangian Analysis of Continua (FLAC) technique for solving the conservation equations for mass and momentum and energy, the heat advection-diffusion equation with frictional heating. For details, see Method section #2. We setup the computational domain configured with dimensions of 1000 km in length and 300 km in depth, comprising 550 by 113 quadrilateral elements. The finest grid spacing, situated at the top center of the domain where the trench is located, is 1 km by 0.5 km. The initial model set up is shown in Fig. S2 and the parameters for initial model are shown in Table S3. The tectonic units, arranged from left to right, include the Eurasian continental shelf, continental slope, continent-ocean transition, South China Sea (SCS), and Philippine Sea Plate. The initial crustal thickness for each block is referenced from seismic reflection images at the northern SCS margin and crossing the Hengchun Peninsula, where deformation is minimal, making these regions suitable proxies for the initial condition. A short slab of oceanic crust, approximately 6 km thick and 70 km long, is introduced into the mantle as a mechanical heterogeneity to facilitate subduction initiation. This configuration is a shortcut to achieve stable subduction quickly. A more realistic subduction initialization model would be proceeded by a slow and long convergence. For boundary velocity, the Luzon Arc advances towards the Eurasian continent at a rate of 7-8 cm/yr along a 300°-310° azimuth. The effective E-W convergence rate at the Manila Trench is approximately 5 cm/yr, which is equally partitioned between the two plates (i.e., 2.5 cm/yr each side) to maintain the trench at the center of the model and within the high-resolution region. Various convergence partitioning schemes tested in uniform-resolution models indicate that the partitioning does not notably impact the model outcomes. The bottom boundary of the model is implemented as a Winkler foundation, allowing for both inflow and outflow, while the top boundary is a free surface incorporating hillslope diffusion to simulate surface process-erosion and sedimentation-controlled by a diffusivity coefficient (κ), which we set as here. To consider the thermal structure for the geodynamic models in the seismic source area in the northern Manila subduction zone, we follow the temperature in the literature. Magnetic anomalies have been used to derive Curie temperature depths ( ~ 550 °C), suggesting that the Curie depth lies approximately 10 km beneath the Moho in the western offshore region of Taiwan. A thermal structure model based on steady-state heat conduction, constrained by near-surface measurements of thermal conductivity and radiogenic heat production, further characterizes the regional thermal regime. This model shows that the Moho temperature in the southern western offshore of Taiwan is about 500 °C. Additionally, a thermal model for the northern Manila Trench, constrained by heat flow observations, indicates that the temperature at the Moho is approximately 150-250 °C in the Philippine Sea Plate and 450-550 °C on the subducting thinned continental slab. We explore several potential configurations of lithospheric rheology, which is not directly observable, and evaluate their predicted strain rates, stress regimes, accumulated plastic strain, and morphology to identify the model that best aligns with available observations. These observations include subduction zone morphology, Moho temperature, lithosphere-asthenosphere boundary (LAB) depth, crustal thickness, constrained by seismic tomographic images, as well as shear strain distributions inferred from seismic azimuth anisotropy.