The Pacific plate, which is the largest oceanic plate on Earth, has implications for the general understanding of plate dynamics, including the origin of intraplate stress and the driving force for plate motion. However, this is currently limited by the scarcity of geophysical and geological observational data. In this study, an instantaneous global mantle flow calculation was performed to predict the intraplate stress field and stress regimes on the Pacific plate using a geodynamic model based on the density anomaly structure of the mantle converted from a seismic tomography model incorporating subducting plates. The numerical results demonstrate that the southern part of the Pacific plate is dominated by a normal faulting regime. In contrast, the northern part is dominated by a thrust faulting regime, in which the tensional stress axes in the older and stable part of the Pacific plate tend to be oblique to the direction of plate motion. This suggests that the stress state of the Pacific plate is almost neutral (i.e., neither compressional nor tensional) along the direction of plate motion. Furthermore, shallow positive buoyancy-induced asthenospheric flow is essential for reproducing the observed plate motion of the Pacific plate.

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 [24] 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 [714], 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 [1719]. 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.

Following previous work (e.g., [20]), the dimensionless conservation equations for mass and momentum with an infinite Prandtl number under the Boussinesq approximation were as follows:
(1)v=0,(2)p+ηv+vTr+Raiζ3δρer=0,
where v is the velocity, p is the dynamic pressure, η is the viscosity, δρ is the density anomaly, er is the unit vector in the radial direction, and Tr indicates the tensor transpose. The dimensionless parameters Rai and ζ are the “instantaneous Rayleigh number” [15] and the mantle/shell radius ratio, respectively, which are defined as follows:
(3)Raiρ0gb3η0κ0,ζbr1,
where ρ0 is the reference density, g is the gravitational acceleration, η0 is the reference viscosity, κ0 is the reference thermal diffusivity, r1 is the Earth’s radius, and b is the thickness of the mantle (note that Rai for this instantaneous mantle flow calculation is not the usual Rayleigh number used in time-dependent mantle convection models). When ρ0=3,300kgm3, g=9.81ms2, η0=1021Pas, κ0=106m2s1, r1=6,371km, and b=2,891km, which was estimated as Rai~7.822×108 and ζ~ 0.454.
The density anomaly (δρ) and temperature anomaly (δT) in the mantle were converted from the S-wave velocity anomaly (δVs) structure from the “TX2019slab” seismic tomography [21] (hereafter, I call it the “TX2019” model, because slab configuration model incorporated in the TX2019slab was not used, as stated below). These conversion equations are expressed as follows:
(4)δρ=αdASdδVS=δρ/TδTδVS/TδTδVS=BSdδVS,(5)δT=δVST1δVS=AS1dδVS,
where ASd is the temperature derivative of the S-wave seismic velocity anomaly as a function of depth (Figure 2(a)), αd is the coefficient of thermal expansion (Figure 2(b)), and BSd=αd/ASd is the conversion factor (Figure 2(c)). These profiles were obtained from the results of mineral physics studies [22]. The “ak135-F” model was used as the reference model of S-wave velocity [23], the temperature of the bottom of the lithosphere (i.e., the temperature of the asthenosphere), T1, was set at 1,330°C [24], and the superadiabatic temperature difference across the mantle was set at 2,000°C [25].

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α0T1T0/2 [29], where T0 is the temperature at the top surface, 0°C. Thus, Δρslab=65.8kgm3. 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).

The viscosity of the mantle varied vertically and laterally (Figures 3(d)–3(f)). In the model, the mantle was composed of four layers—highly viscous lithosphere, low-viscosity asthenosphere with a bottom depth of 220 km, upper mantle (with depths of 220–660 km), and lower mantle (with depths of 660–2,891 km). Although previous mineral physics (e.g., [30, 31]) and geodynamic modeling (e.g., [32, 33]) studies have suggested a depth-dependent viscosity for the entire lower mantle, it remains a controversial issue. Thus, the present numerical model assumes a single layer of the lower mantle for simplicity [15, 3438]. The reference viscosity of the upper mantle was set as η0=1021Pas. In addition, the viscosity of the mantle (except for the lithosphere) depends on the temperature T=T1+δT, as follows:
(6)ηT,Γ=ΔηΓexp2T1EaRT+T1Ea2RT1,
where R is the gas constant (8.3144 J K−1 mol−1), Ea is the activation energy (335 kJ mol−1) [39, 40], and ΔηΓ is the viscosity contrast between each layer Γ and the upper mantle. The viscosity contrast between the lower and upper mantle viscosities was set to ΔηΓ=Δηlwm=30 for the reference model and 100 in some other models (e.g., [38]). The viscosity defined in Equation (6), namely, the viscosity of temperature anomaly converted by Equation (5) is truncated from 1018 Pas to 1024 Pa s for the numerical stability.
The thickness of the oceanic lithosphere was given by [25]
(7)dlit=2.32κ0ap,
where κ0 is the thermal diffusivity (10−6 m2 s−1), ap is plate age given by a global crustal age model [41], and the maximum bottom depth of the oceanic lithosphere was set as dlitmax=100km. The spatial distribution of the continental lithosphere was given by the “TC1” model, i.e., the thermal model of the continental lithosphere [42]. The viscosity of the lithosphere (ηlit) is not followed by Equation (6) and, instead, logarithmically decreased with depth to be η0=1021 at the bottom of the lithosphere:
(8)log10ηlitd=log10ηlitmaxlog10ηlitmaxlog10η0dlitmaxd,
where d is the depth. The maximum viscosity at the top surface boundary ηlitmax is fixed at 3×1023Pas for the reference model (e.g., [43]) and 3×10221023Pas for some other models. In the present study, the viscosity contrast between the lithosphere and upper mantle, ΔηΓ=Δηlit=ηlitmax/η0, and that between the asthenosphere and the upper mantle, ΔηΓ=Δηast, were taken as free parameters in this model (see Section 3). The viscosity of the subducting slab was assumed to equal Δηlit. The density anomaly of the craton is assumed to be zero for simplicity, whereas the viscosity of the craton was assumed to be one order of magnitude larger than that of the surrounding lithosphere, as suggested by a geodynamic model featuring craton longevity [43].
Furthermore, mechanically weak plate margins (or weak plate boundaries) were imposed in the lithosphere from the data used in the Global Strain Rate Map (GSRM) version 2.1 [44]. The viscosity of the plate boundary was estimated as
(9)ηpbθ,ϕ=τIIε̇IIθ,ϕ,
where ε̇II is the second invariant of deviatoric strain rate at the surface from the GSRM dataset, and the deviatoric stress was fixed at τII=3MPa [45]. Thus, the obtained viscosity of the plate boundaries depends on the spatial distribution of ε̇IIθ,ϕ (see Figure 2 in [44]). The lateral viscosity variations of the plate boundaries were determined by Equation (9), whereas the depth-dependent viscosity of the lithosphere is determined by Equation (8) (Figure 3(d)). The plate boundaries penetrated into all the depths of the lithosphere.
After the velocity and stress fields were obtained by the instantaneous mantle flow calculation by solving conservation equations for mass and momentum (Equations (1) and (2)), the maximum and minimum horizontal principal deviatoric stress axes and the vertical principal deviatoric stress axis (i.e., σHmax, σhmin, and σv, respectively) were analyzed in each computational grid. First, the three faulting regimes, that is, normal faulting (NF), strike-slip faulting (SS), and thrust faulting (TF) regimes, were defined by the magnitudes of σHmax, σhmin, and σv, respectively. Following the definition from previous work (e.g., [1, 46]) and the World Stress Map project (e.g., [2]), the regions of NF, SS, and TF were classified into three “stress regimes” that satisfy the following conditions:
(10)σv>σHmax>σhminforNFregime,σHmax>σv>σhminforSSregime,σHmax>σhmin>σvforTFregime,
where σHmax and σhmin are the maximum and minimum horizontal principal deviatoric stresses, respectively, and σv is the vertical principal deviatoric stress (Table 1).
Then, each of the three faulting regimes was further classified into three stress regimes based on the magnitude of the stress ratio, Rs, derived as follows:
(11)Rs=σ2σ3σ1σ3,
where σ1, σ2, and σ3 are the maximum (i.e., compressional), medium, and minimum (i.e., extensional) principal deviatoric stresses, respectively. Any of these stresses corresponds to σHmax, σhmin, or σv (e.g., [47]). Similar to the work of Fitzenz and Miller [48], the three faulting regimes were further classified into three stress regimes depending on the stress ratio Rs, which is continuously ranged from 0 to 1. Following previous work [49], the criteria were set at 0Rs1/3, 1/3Rs2/3, and 2/3Rs1 (Table 1). For the normal faulting regime, the portion of the strike-slip component increases with increasing Rs, whereas, for the thrust faulting regime, the portion of the strike-slip component increases with decreasing Rs. For the strike-slip faulting regime, the portion of the thrust faulting component increases with decreasing Rs, and the portion of the normal faulting component increases with increasing Rs (see also Yoshida [49] for details).

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.

3.1. Reference Model

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−1 and 4.2 cm yr−1 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−1) [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−1, 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−1 (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.

3.2. Validation of the Predicted Plate Motion

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.

The relative error is defined in terms of global root-mean-square (RMS) misfit as follows (e.g., [57]):
(12)Err=vpvo2d2Svo2d2S,
in which vp and vo are the vector fields of the predicted and observed plate motions, respectively, and d2S represents integration over the target surface area, S. Based on Err, the variance reduction is defined by
(13)VR=1Err2.

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.

3.3. Effects of Shallow Positive Buoyancy-Induced Asthenospheric Flow

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.

3.4. Toroidal-Poloidal Ratio

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.

4.1. Major Driving Forces of the Pacific Plate

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.1cmyr1 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.3cmyr1 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, 7275]). 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.

4.2. Effects of Low-Density Anomalies in the Shallow Mantle

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.

4.3. Slab Pull Force and Numerical Limitations

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.

Appendix

Toroidal-Poloidal Ratio

The ratio of the power spectra of the toroidal and poloidal motions (i.e., the T/P ratio) at each depth for each model was expressed as follows (e.g., [49, 66, 80]):
(A.1)ΓT/Pr=llmaxm=0latorlmc2+atorlms2llmaxm=0lapollmc2+apollms2×100%,
where l and m are the degree and order of the spherical harmonic function, lmax is set to 16, and aporimi and atorimi (where superscript i is the sine or cosine component) are the coefficients of the poloidal and toroidal components, respectively.
(A.2)apollmir=14πll+1SvθYlmθ+vϕsinθYlmϕdS,atorlmir=14πll+1SvθsinθYlmϕvϕYlmθdS.
Here, Ylm is the fully normalized spherical harmonic function (e.g., [81]), vθ and vϕ are the latitudinal and longitudinal velocities, respectively, and S is the entire surface area at each depth. The T/P ratio at the top surface, γT/P=ΓT/Pr1, for each model is listed in Table 2. Note that Equation (A.1) at the top surface (that is defined in the spherical harmonic spectral domain) represents the spatial measure of relative toroidal/poloidal RMS velocity amplitudes as follows:
(A.3)T/P=vtor2d2Svpol2d2S×100%,
where vtor and vpol are the toroidal and poloidal components of the surface velocity field, respectively.

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

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).