Stress State of the Stable Part of the Pacific Plate Predicted by a Numerical Model of Global Mantle Flow Coupled with Plate Motion

The intraplate stresses of the Earth’s oceanic plates, which can be measured by geophysical and geological approaches based on focal mechanisms, borehole breakouts, in situ stress measurements, fault slip analysis, and volcanic vent alignments, are crucial for understanding the origin of the intraplate stress field and the driving force of plate motion (e.g., [1]). The Pacific plate, which is the largest oceanic plate on Earth, has implications for achieving a general understanding of such plate dynamics. However, because there are a few large seismicities and thus few observational datasets for the stable part of the Pacific plate far from the spreading ridges and trenches, it remains difficult to infer the intraplate stress characteristics. Indeed, reliable, high- to moderate-quality data of focal mechanisms and breakouts in the stable part of the Pacific plate listed in the World Stress Map (WSM) database [2–4] are sparse compared to the large size of this plate.

Figure 1 demonstrates the distribution of the faulting regimes in the Pacific plate from the WSM database. The directions of bars for the normal, strike-slip, and thrust faulting regimes indicate the orientations of maximum horizontal stress axes (⁠$σHmax$⁠). Based on the limited focal mechanism data for the Pacific plate, normal faults are dominant near the East Pacific Rise (EPR), whereas thrust and transform faults are dominant in the older parts far from the EPR [5, 6]. However, the “quality” of almost all of the data in the stable part of the Pacific plate is classified as “C,” i.e., moderate quality, not higher qualities (“A” and “B”) (Figures S1 and S2 in the Supplementary Materials). The horizontal maximum principal stress axes with parallel and perpendicular orientations along the spreading ridge under the thrust faulting regimes seem to have coexisted near a part of the spreading ridge (black circles marked “A” and “B” in Figure 1 and see the close-up view in Figure S1 in the Supplementary Material). On the other hand, there is no tendency for the orientations of $σHmax$ in the stable part of the Pacific plate, which presumably implies that these data do not necessarily represent the stress field of the Pacific plate derived from global mantle flow circulation, but rather from local tectonics. Furthermore, few relevant data have been obtained from the South Pacific region, in which hotspots are concentrated along some hotspot chains (the region enclosed in “C” and “D” in Figure 1).

Assuming that the direction of $σHmax$ tends to be parallel to the direction of plate motion, the Pacific plate far from the EPR is generally considered to be dominated by compressional stress along the direction of plate motion. A previous systematic numerical study of instantaneous global mantle flow demonstrated that from a best-fitting model that well explains the plate velocity and strain rate of Earth, the Pacific plate near the EPR is dominated by a normal faulting regime with compressive horizontal principal deviatoric stress axes almost perpendicular to the direction of plate motion, whereas the older part of the Pacific plate is dominated by the strike-slip regime with compressive horizontal principal deviatoric stress almost perpendicular to the direction of plate motion [7]. These results indicate that the stable part of the Pacific plate is subject to tensional stress along the direction of plate motion.

Although the global stress field and driving force of plates have been a subject of interest in geodynamics and studied using semianalytical calculations of instantaneous global mantle flow [7–14], it has not been focused on the Pacific plate with few observation data. This study addresses the stress state of the entire Pacific plate based on an instantaneous global mantle flow calculation with a full formulation of the momentum equation used in numerical simulations of time-dependent mantle convection. The numerical model is based on the density anomaly structure of the entire mantle converted from a recent high-resolution seismic tomography model incorporating subducting slabs. The intraplate stress field of the Pacific plates is quantitatively evaluated in terms of the stress ratio and “stress regimes” using velocity and stress fields obtained from the global mantle flow model with plate motion and subducting plate.

A mantle flow model confined in a three-dimensional (3D) spherical shell domain along spherical polar coordinates (⁠$r$⁠, $θ$⁠, $ϕ$⁠) was used to obtain the intraplate stress field in the global domain. The mantle, with a thickness of 2,891 km, was modeled as a viscous fluid with an infinite Prandtl number, and the instantaneous mantle flow was numerically calculated using the mantle convection code “ConvGS” based on the finite-volume method (e.g., [15, 16]) and the Yin-Yang grid [17–19]. The number of finite volumes was $128r×128θ×384ϕ×2$⁠. The mechanical conditions for the top and bottom surface boundaries of the mantle were impermeable and shear stress-free.

Similar to previous work [26], the density anomaly in the regions with high-seismic-velocity anomalies above a depth of 300 km—except for the regions of subducting slab—was set to zero (i.e., $δρθ,ϕ=0$⁠) to remove compositional effects in the shallower part of the upper mantle. Instead, the configuration of subducting slabs was constructed from a global slab configuration model, “Slab2” [27], which includes a higher resolution of slab geometry in each subduction zone compared to a previous slab model, “Slab1.0” [28]. The density anomaly of the subducting slab was given by $Δρslab=ρ0α0T1–T0/2$ [29], where $T0$ is the temperature at the top surface, 0°C. Thus, $Δρslab=65.8kgm−3$⁠. The density anomaly structures at depths of 373 km and 2,800 km are shown in Figures 3(a) and 3(b), respectively. The 3D isosurface image of the density anomaly of the whole mantle is shown in Figure 3(c).

Figure S3 in the Supplementary Materials shows the distribution of the horizontal principal axes of the deviatoric stress for the entire surface in geometrically symmetric, steady-state mantle convection with a low Rayleigh number [49]. This indicates that a normal faulting regime is dominant on six cylindrical upwelling plumes, whereas a thrust faulting regime is dominant on sheet-like downwelling plumes, as is the case for real subducting plates. It is noteworthy that in this paper, the figures showing the distribution of the horizontal principal axes of deviatoric stress, the two horizontal axes (i.e., $σHmax$ and $σhmin$⁠) are displayed while $σHmax$ in the pure normal faulting regime (shown by “N”) and $σhmin$ in the pure thrust faulting regime (shown by “T”) are omitted to visually distinguish these regimes.

The studied models are listed in Table 2. The numerically obtained plate motion (i.e., “predicted” plate motion) was compared with the observed plate motion of the NNR-MORVEL56 model in the no-net-rotation reference frame [50] (Figure 4). The viscosity contrast between the lithosphere and upper mantle was set to $Δηlit=300$ for the reference model (Model Ref) based on a previous numerical study [51] and the preliminary runs of the present study. When $Δηlit=300$⁠, the horizontal root-mean-square (RMS) velocity at the top surface, $Vh_RMS$⁠, was 5.1 cm yr⁻¹ and 4.2 cm yr⁻¹ for Model Ref (⁠$Δηlwm=30$⁠) and Model Ref_Lwm100 (⁠$Δηlwm=100$⁠) (Table 2), respectively, which is comparable to the RMS velocity of the observed plate motion (4.4 cm yr⁻¹) [50]. The surface RMS velocity increased when $Δηlit$ decreased, as expected in Figure S4A in the Supplementary Materials). Incidentally, when the density anomalies in the mantle were converted from P-wave velocity anomaly structures from the TX2019 model (Model Ref_P), the speed of the Pacific plate was too slow and the $Vh_RMS$ was 2.3 cm yr⁻¹, which is approximately half of the observed plate motion RMS velocity (Figure S5 in the Supplementary Materials). This is due to the smaller magnitude of the seismic velocity perturbations in the P-wave seismic tomography model that results in the smaller magnitude of the density anomaly variations to drive the mantle flow.

The viscosity contrast between the asthenosphere and the upper mantle, $Δηast$⁠, was set to 0.01 for the reference model and increased by up to 1 for the other models. Based on the preliminary analyses, the horizontal RMS velocity in the entire top surface of the models with $Δηast=0.1$ was low compared with the models with $Δηast=0.01$ (Table 2), and the pattern of the intraplate stress in the Pacific plate did not significantly depend on the magnitude of $Δηast$⁠. When $Δηast=1$⁠, i.e., there is no low-viscosity asthenosphere, the Pacific plate is driven by deeper mantle flow below the asthenosphere, and the horizontal RMS velocity at the top surface decreased by 1.6 cm yr⁻¹ (Table 2), comparable to the underlying speed of the mantle flow (Figure S4C in the Supplementary Materials).

Figure 5 shows the numerical results for the reference model (Model Ref) with a viscosity contrast between the lithosphere and the upper mantle of $Δηlit=300$ and a viscosity contrast between the lower and upper mantle of $Δηlwm=30$⁠. The three panels are included in Figure 5, i.e., the distribution of the viscosity and horizontal velocity fields near the entire top surface boundary (Figure 5(a)), the distribution of the stress ratio $Rs$ in each faulting regime near the entire top surface boundary (Figure 5(b)), and the distribution of the horizontal principal deviatoric stress axes in the Pacific plate (Figure 5(c)). It was found that the direction and speed of the Pacific plate are comparable to the observed plate motions (Figure 5(a)). The distributions of $Rs$ (Figure 5(b)) and the horizontal principal deviatoric stress axes (Figure 5(c)) demonstrate that normal faulting and strike-slip faulting regimes are broadly dominant in the southern half of the Pacific plate. In particular, in the French Polynesia region (Figure 1(D)), where some hotspots with large buoyancy fluxes have been observed (yellow triangles in Figure 5(c)), the compressional stress is generally dominant in the east-west direction. In contrast, the tensional stress is dominant in the north-south direction (Figure 5(d)), which is consistent with previous result [7].

In contrast, the northern part of the Pacific plate is dominated by a thrust faulting regime in which the direction of tensional stress axes in the stable and older parts of the Pacific plate tends to be oblique along the direction of plate motion (Figure 5(e)). This result implies that the stress state is nearly neutral (i.e., neither compressional nor tensional) along the direction of the plate motion. On the other hand, the compressive horizontal principal deviatoric stress axes are nearly perpendicular to the Japan and Izu-Ogasawara Trenches (Figure 5(f)), which is broadly consistent with the previous results [7].

The present numerical studies suggest that the entire Pacific plate inherently maintains a neutral stress state as shown in both the northern and southern parts of the plate. However, there is a dichotomy of stress regimes between the two hemispheres: the thrust faulting regime is dominant in the northern Pacific, whereas the normal faulting regime is dominant in the southern Pacific (Figure 5(b)). This may be attributed to large-scale upwelling flow originated from the “Large Low Velocity Provinces (LLVPs)” [52] under the South Pacific, as shown in the density anomaly model converted from seismic velocity anomalies (Figures 3(b)–3(c)).

The background gray scales and the length and width of the arrows shown in Figure 5(c) demonstrate that the square root of the second invariant of deviatoric stress, $τII$⁠, is generally in the order of 10 MPa and not larger than 100 MPa in the stable part of the Pacific plate, whereas the viscosity of weak plate boundaries is lower than the order of 1 MPa. These stress magnitudes are comparable to those in the stable part of the plates and plate boundaries of the Earth [53] and are not significantly different from those even for the other models studied in this study.

In this study, the “relative error” (⁠$Err$⁠) and “variance reduction” (⁠$VR$⁠) were used to evaluate the fit of the predicted velocity field to the observed plate motions. Because intraplate stress is a focus of this study, the plateness—which is required to determine the Euler pole of single plate motion based on the assumption of rigidity—was not calculated as in previous studies (e.g., [54, 55, 56]). However, it is sufficient to only assess the relative misfit and relevance among the models studied in this paper.

In this study, $Err$ and $VR$ are analyzed over (i) the entire surface of the Earth, (ii) the entire surface of the Pacific plate, and (iii) the entire region outside of the Pacific plate (Table 2).

In the Model Ref (Figure 5), $Err$ and $VR$ over the entire Earth’s surface are 0.57 and 0.68 (Table 2), respectively, which indicates that the prediction of the observed plate motion is relatively valid for the finite-volume approach used in this study, compared with geodynamic models based on the torque balance method and finite-element approach (see discussion in Section 4.3). When confined to the Pacific plate, $Err$ and $VR$ are 0.14 and 0.98, which indicates better validity compared to the method over the entire surface. When the viscosity of the lower mantle is increased to $Δηlwm=100$ for Model Ref_Lwm100 (Figure S6 in the Supplementary Material), the $Err$ and $VR$ are improved to 0.50 and 0.75 over the entire Earth’s surface and 0.10 and 0.99 in the Pacific plate, respectively (Table 2), which indicates that Model Ref_Lwm100 is the “best” among the models studied in this study. Figure 6 shows that because of the low viscosity in the asthenosphere, the flow speed in the top part of the asthenosphere is larger than that of the plate motions for both Models Ref (Figure 6(a)) and Ref_Lwm100 (Figure 6(b)). Thus, it is important to study the mechanical coupling between the plate motion and asthenospheric flow (see Section 3.3) and the possible driving forces behind the plate motion (see Section 4.2).

In models with contrasting viscosities between the lithosphere and the upper mantle of $Δηlit=100$ (Model Ref_Lit100) and $Δηlit=30$ (Model Ref_Lit030), dichotomous stress regimes were seen between the two hemispheres throughout the entire Pacific plate (see Figures S7 and S8 in the Supplementary Material). However, the validity in terms of the prediction of the observed plate motion became too low (i.e., $VR$ over the Earth’s surface and the Pacific plate were below 0.23 and 0.47, respectively; see Table 2), because the speed of the Pacific plate becomes too large compared with that of the model with $Δηlit=300$⁠.

Previous geodynamic models have suggested that asthenospheric flow induced by the low-density anomaly beneath the lithosphere plays an important role in the formation of laterally moving large-scale plates (e.g., [58, 59]), and mantle drag by the underlying asthenospheric flow can act as one of the primary driving forces for plate motion under limited conditions (e.g., [60]). To observe the effects of shallow positive buoyancy-induced asthenospheric flow on the motion and intraplate stress of the Pacific plate, additional models were tested, namely, Model NLD where $Δηlwm=30$ and Model NLD_Lwm100 where $Δηlwm=100$⁠. In these models, the low-density anomaly regions converted from the low seismic velocity anomaly above a depth of 300 km (Figure 7(a)) were replaced by a no-density anomaly and no-temperature anomaly, i.e., $δρθ,ϕ=δTθ,ϕ=0$ (Figure 7(c)), compared with the model which considers density anomalies in the uppermost upper mantle on a global scale (e.g., [61]) and in a regional model (e.g., [62, 63]). The lateral viscosity variations owing to the dependence of temperature anomaly (Equation (6)) were also eliminated above a depth of 300 km (cf. Figures 7(b) and 7(d)).

Figure 8 demonstrates that when there is no asthenospheric flow driven by shallow low-density anomaly, the motion speed of the Pacific plate decreases (Figure 8(a)) relative to that of Models Ref and Ref_Lwm100. This is because the horizontal speed of the asthenospheric flow that drives the plate motion was reduced (Figure 9(a)), due to the removal of density anomalies in the asthenosphere, but also the asthenosphere becoming more viscous for both Models NLD (Figure 6(c)) and NLD_Lwm100 (Figure 6(d)), compared with Models Ref (Figure 6(a)) and Ref_Lwm100 (Figure 6(b)). This implies that active asthenospheric flow originating from the low-density anomaly in the shallow part of the mantle is important for reproducing the observed plate motion of the Pacific plate. In Models NLD and NLD_Lwm100 (Figure S9 in the Supplementary Material), the $VR$ over the entire Earth’s surface was 0.57 and 0.44 and 0.62 and 0.47 in the Pacific plate, respectively, which indicates a lower validity of the prediction of the observed plate motion (Table 2).

Remarkable differences between the results of Model Ref and Model NLD were observed in the patterns of the stress ratios and stress regimes on the Pacific plate. Specifically, in Model NLD, the distributions of the stress ratio and stress regime show that a normal faulting regime is broadly dominant in the southern part of the Pacific plate, whereas a strike-slip faulting regime is broadly dominant in the northern part of the Pacific plate (Figure 8(b)). On the other hand, for the northern part, the principal deviatoric compressional and tensional stress axes tend to be oblique along the direction of plate motion (Figures 8(c) and 8(d)), which is similar to the results of Model Ref.

The ratio of the power spectra of the toroidal and poloidal motions (i.e., the T/P ratio) in the lithosphere is a measure of the plate-like feature on the modeled Earth (see the appendix). It is accepted that the toroidal component of the observed plate motion is comparable to the poloidal component (i.e., “toroidal-poloidal equipartitioning”) [64, 65]. Namely, the T/P ratio of the observed plate motion at the surface is ~100%, which is supported by the previous instantaneous global mantle flow model, in which the viscosity of the plate boundaries is relatively weak compared to the stable portion of the plate [66]. Furthermore, a more detailed analysis of the T/P ratio as a function of each spherical harmonic degree using the observed plate motion data from the NNR-MORVEL56 model supports the hypothesis of toroidal-poloidal equipartitioning [67].

The T/P ratio at the top surface, $γT/P=ΓT/Pr1$⁠, for each model is listed in Table 2, which demonstrates that $γT/P$ is in the range of approximately 31 to 55% in the models of the present study, and as far as Models Ref, Ref_Lwm100, NLD, and NLD_Lwm100, γ_T/P ranges from ~43% to 55%. This range is consistent with a previous systematic numerical model of time-dependent mantle convection, with self-consistent generation of the surface plate motion in 3D Cartesian geometry [68] and 3D spherical-shell geometry [69, 70]. Meanwhile, this $γT/P$ range is lower than the observed T/P ratio, likely because the average viscosity of the plate boundaries in the present model is higher than that of the actual plate boundaries. Yoshida et al. [66] performed a systematic numerical study and determined that the T/P ratio becomes high (close to ~100%) as the viscosity ratio between the stable part of the plates and the plate boundaries increases. However, the relative difference in the T/P ratio between the models is useful for discussing the magnitude of the absolute toroidal energy in each model.

Figure 9 demonstrates the RMS horizontal velocity at each depth (Figure 9(a)) and the ratio of toroidal and poloidal energy at each depth (Figure 9(b)) for Models Ref, Ref_Lwm100, NLD, and NLD_Lwm100. It is found that the laterally averaged speed of the asthenospheric flow was decreased, whereas the T/P ratio increased when shallow positive buoyancy-induced asthenospheric flow was removed (Figure 9(b)). This is because the toroidal energy originating from the lateral viscosity variations in the mantle is not well released in the asthenosphere and the upper mantle but, instead, is released in the lithosphere (see also the small subpanel in Figure 9(b)). For the two studied models (i.e., Models Ref and NLD), increasing the viscosity of the lower mantle, i.e., $Δηlwm=100$ (Models Ref_Lwm100 and NLD_Lwm100), did not significantly affect the patterns of plate motion and intraplate stress (Figures S6 and S9 in the Supplementary Materials) nor the T/P ratios in the asthenosphere and lithosphere (Figure 9(b)). This implies that the dynamics of the shallower part of the mantle are not affected by lower mantle viscosity, although the laterally averaged speed of the asthenospheric flow was slightly decreased relative to Model Ref with $Δηlwm=30$ (Figure 9(a)).

Figure 10 shows the numerical results for Models Ref_Ast0.1 with $Δηast=0.1$⁠. The viscosity contrast between the lithosphere and the upper mantle and the viscosity contrast between the lower and upper mantle are the same as Model Ref (Figure 5), i.e., $Δηlit=300$ and $Δηlwm=30$⁠. Increasing the viscosity of the asthenosphere, i.e., $Δηast=0.1$ (Models Ref_Ast0.1 and NLD_Ast0.1) or $Δηast=1.0$ (Models Ref_Ast1.0 and NLD_Ast1.0), respectively, also did not significantly affect the patterns of plate motion and intraplate stress (cf. Figures 10 and S10 in the Supplementary Materials and Figures S11 and S12 in the Supplementary Materials). The T/P ratio in the lithosphere varied by approximately 10% among the models with different $Δηast$ (Table 2 and Figure S4D in the Supplementary Materials) and $Δηlit$ (Table 2 and Figure S4B in the Supplementary Materials). As expected, the $γT/P$ increased with a decrease of the asthenospheric viscosity and, thus, increase of the lateral flow velocity in the asthenosphere; the $γT/P$ values are 43.1%, 35.2%, and 31.3% for Models Ref with $Δηast=0.01$⁠, Models Ref_Ast0.1 with $Δηast=0.1$⁠, and Models Ref_Ast1.0 with $Δηast=1.0$⁠, respectively. On the other hand, the $VR$ of the entire Earth’s surface in Model Ref_Ast0.1 with $Δηast=0.1$ is 0.70, which is a relatively higher value than those for Model Ref with $Δηast=0.01$ (⁠$VR=0.68$⁠) and Model Ref_Ast1.0 with $Δηast=1.0$ (⁠$VR=0.47$⁠) (Table 2). This result may be consistent with the findings of a previous report, which states that the reasonable viscosity difference between the asthenosphere and the upper mantle is around 0.1, as suggested by a previous geodynamic model of craton longevity [43]. Namely, if the viscosity of the asthenosphere is too low or close to 0.01, the cratonic lithosphere becomes mechanically weak and cannot ensure longevity over Earth’s history.

The instantaneous global mantle flow calculation in the present study demonstrated that shallow positive buoyancy-induced asthenospheric flow is essential for reproducing the observed plate motion of the Pacific plate. To further address this conclusion, an additional model was studied, in which low-density anomalies above a depth of $dNLD=200km$ and 100 km were removed (Models NLD_Dep200 and NLD_Dep100), compared to Model NLD which removed these anomalies at 300 km. In the model with $dNLD=200km$⁠, the speed of the Pacific plate becomes quite slow (⁠$VRMS=2.1cmyr−1$ and $VR=0.66$ in the Pacific plate) compared with the observed plate motion (Figure S13 in the Supplementary Materials). In contrast, when $dNLD$ is changed to 100 km, the predicted plate motion is somewhat improved (Figure S14 in the Supplementary Materials), as revealed with $VRMS=4.3cmyr−1$ and $VR=0.99$ in the Pacific plate. These findings, along with the results for Model Ref (Figure 5), imply that the active shallow positive buoyancy-induced asthenospheric flow beneath the plate is crucial for fully characterizing the entire patterns of plate motion and intraplate stress. Several authors have suggested that for all plates on a global scale, the ridge push force is only a minor contribution to the surface plate motions (e.g., [71, 72–75]). On the other hand, for the Pacific plate, our numerical results showed that shallow positive buoyancy-induced asthenospheric flow is essential for reproducing the observed plate motion of the Pacific plate.

Several previous papers have confirmed the dynamic importance of deep-seated positive buoyancy (>250 km) under the EPR as a major driving force of Pacific plate motion through the asthenospheric drag. For example, Forte et al. [76] found that the slowing of the Pacific and the Nazca plates is caused by the time dependence of the buoyancy of the western Pacific slab and the nonsteady decrease in the deep positive buoyancy below the EPR. Their mantle convection model successfully reproduced current plate velocities and retained the gravity and topography of the Earth’s surface. Glišović and Forte [57] concluded that based on a reconstruction of mantle evolution during the Cenozoic, an active upwelling associated with a dynamically stable mantle plume from the core-mantle boundary has been maintained below the EPR for 65 million years. This finding supports the hypothesis that the EPR is directly forced by a long-lived and stable mantle plume and explains the observation that the EPR ridge has been laterally stable over the past 80 million years. Furthermore, Rowley et al. [77] concluded that deep-seated mantle upwelling beneath the EPR drives the horizontal component of the asthenosphere beneath the plate, which is faster than the surface plate, which contributes to the motion of the Pacific plate through the asthenospheric drag.

Note that the present study focused on the effects of removing low-density anomaly above the depth of 300 km on the results; this is different from Rowley et al. [77] that focused on the deep-seated buoyancy at the depths of >250 km or even the lower mantle below EPR. Rowley et al. [77] demonstrated that removing the deep-seated positive buoyancy below the EPR primarily affected the magnitudes instead of the directions of the Pacific plate. On the other hand, the analysis using the present model showed that removing the shallow positive buoyancy (<250 km) also considerably affects the direction of the Pacific plate motion. However, the present study removed the shallow low-density anomalies globally instead of only below the EPR, which means that the shallow positive buoyancy below the EPR is only part of global shallow positive buoyancy, and thus, it does not adequately distinguish the effects of low-density anomaly below the EPR and that in a global scale. Nevertheless, Figures 6(a) and 6(b) show that the flow speed originating below the EPR seems to be large in the entire domain, implying that the effects of low-density anomalies below the EPR are dominant.

In the discussion of Section 4.1, it is difficult to determine whether the modeled differences in plate motion are due to the removal of the low-density anomaly in the shallow part of the mantle or the removal of the lateral viscosity variation based on Equation (6). Here, to assess the effects of the low-density anomaly in the shallow mantle alone on the change in the plate motion of the Pacific plate, four additional models were studied: Models Ref_NLV, NLD_YLV (cf. Figures 7(e) and 7(f)), NLD_Dep200_YLV, and NLD_Dep100_YLV (Figures S15, S16, S17, and S18 in the Supplementary Materials, respectively).

When comparing the reference model (Model Ref) with Model Ref_NLV, in which low-density anomalies at depths above 300 km are considered but the effect of lateral viscosity variations is removed, and Model NLD_YLV, in which low-density anomalies at depths above 300 km are removed but the effect of lateral viscosity variations is considered, the $VR$ over the Pacific plate is 0.98 (Ref), 0.99 (Ref_NLV), and 0.66 (NLD_YLV), respectively (Table 2). These results imply that the low-density anomalies in the shallow mantle alone largely affect the plate motion pattern of the Pacific plate.

Table 2 shows that in Models NLD_dep100, NLD_Dep200, and NLD, all of which removed both the low-density anomalies and lateral viscosity variations, the $VR$ over the Pacific plate significantly decreases from 0.99 (NLD_Dep100) to 0.66 (NLD_Dep200) and 0.62 (NLD). When comparing Models NLD_Dep100_YLV, NLD_Dep200_YLV, and NLD_YLV, all of which retain lateral viscosity variations, it also significantly decreases from 0.99 (NLD_Dep100_YLV) to 0.69 (NLD_Dep200_YLV) and 0.66 (NLD_YLV). These results demonstrate that the existence of lateral viscosity variations in the shallow mantle does not make a large difference in the motion of the Pacific plate; thus, low-density anomalies in the shallower part of the mantle play a more important role in reproducing the observed plate motion of the Pacific plate.

If mantle drag by the underlying asthenospheric flow is one of the primary driving forces for plate subduction, and the slab pull force does not act as the primary driving force, the slab pull force may be offset by the slab resistance force acting on the front edge of the subducting plate [78], likely owing to the highly viscous underlying mantle. However, the mantle flow models constructed in this study cannot fit the observed plate motions over other plates well, including other major continental plates (see $VR$ and $Err$ over the entire Earth surface and the entire region outside of the Pacific plate in Table 2). This is confirmed in each panel showing the pattern of predicted plate motions in Models Ref (Figure 5(a)) and Ref_Lwm100 (Figure S6A), even when the $VR$ is maximized over the Pacific plate. Although this study focused on the Pacific plate, global consistency is another important criterion for evaluating the robustness of global mantle flow models. Thus, the lack of applicability to other plates is a limitation of the mantle flow models assessed here.

The present numerical model based on the finite-volume method provides a model setup with weak plate boundaries in the lithosphere with a maximum thickness of 100 km (Figure 3(d)). In contrast, the horizontal thickness of the subducting slab is comparable to that of the weak plate boundaries owing to the limitation of numerical resolution. Thus, the slab pull force is not incorporated in the model, which is beyond the scope of this study, and the present model only includes the effect of mantle drag. This means that the slab pull force is not included as a stress guide from the strong continuous slab to the surface plate, as it is in the force torque balance (e.g., [74, 79]) or the finite-element methods (e.g., [55]). If the stress guide is not incorporated between the subducting slab and the horizontal part of subducting plate, the plate motions are driven purely by buoyancy forces from lateral density anomalies instead of slab pull forces. Although the issue of exact slab strength under the plate-mantle coupling or decoupling system is beyond the scope of this study because I mainly focused on the effect of the basal tractions from mantle flow, future studies should aim to further develop the horizontal resolution as computation power improves.

As stated in Section 1, very little observational data for the stable part of the Pacific plate far from the spreading ridges and trenches, as well as near the spreading ridges, are available. In the present study, the intraplate stress field of the stable part of the Pacific plate was investigated using an instantaneous global mantle flow model. The conclusions are as follows: (1) The shallow positive buoyancy-induced asthenospheric flow is essential for reproducing the observed plate motion of the Pacific plate. (2) The southern part of the Pacific plate is dominated by a normal faulting regime, whereas the northern part is dominated by a thrust faulting regime. (3) The tensional stress axes in the older and stable part of the Pacific plate tend to be oblique to the direction of plate motion, which suggests that the stress state of the Pacific plate is almost neutral along the direction of plate motion.

All data are available in the manuscript and supplementary material. Raw data can be downloaded from doi:10.5281/zenodo.7326028 and provided on request from the corresponding author.

The authors declare that there is no conflict of interest regarding the publication of this paper.

I thank Thorsten Becker for providing observed plate motion data [82], which was used in this paper and previous paper [83]. The TX2019slab model [21] was taken from the IRIS EMC website (http://ds.iris.edu/ds/products/emc-tx2019slab/). The global slab configuration model, Slab2, was taken from USGS website (doi:10.5066/F7PV6JNV) [84]. Figure 1 is created using CASMO online interface via the WSM project website (http://www.world-stress-map.org/casmo/) [4]. Almost all the figures were produced using Generic Mapping Tools [85, 86]. Numerical simulations were performed on the supercomputer facilities of the DA system and Earth Simulator (ES4) at the Japan Agency for Marine-Earth Science and Technology (JAMSTEC). This study was supported by the Japan Society for the Promotion of Science (JSPS KAKENHI, Grant Numbers JP22K03787 and JP20H00200).