Dynamics of the North American Plate: Large‐Scale Driving Mechanism From Far‐Field Slabs and the Interpretation of Shallow Negative Seismic Anomalies

With a small fraction of marginal subduction zones, the driving mechanism for the North American plate motion is in debate. We construct global mantle flow models simultaneously constrained by geoid and plate motions to investigate the driving forces for the North American plate motion. By comparing the model with only near‐field subducting slabs and that with global subducting slabs, we find that the contribution to the motion of the North American plate from the near‐field Aleutian, central American, and Caribbean slabs is small. In contrast, other far‐field slabs, primarily the major segments around western Pacific subduction margins, provide the dominant large‐scale driving forces for the North American plate motion. The coupling between far‐field slabs and the North American plate suggests a new form of active plate interactions within the global self‐organizing plate tectonic system. We further evaluate the extremely slow seismic velocity anomalies associated with the shallow partial melt around the southwestern North America. Interpreting these negative seismic shear‐velocity anomalies as purely thermal origin generates considerably excessive resistance to the North American plate motion. A significantly reduced velocity‐to‐density scaling for these negative seismic shear‐velocity anomalies must be incorporated into the construction of the buoyancy field to predict the North American plate motion. We also examine the importance of lower mantle buoyancy including the ancient descending Kula‐Farallon plates and the active upwelling below the Pacific margin of the North American plate. Lower mantle buoyancy primarily affects the amplitudes, as opposed to the patterns of both North American and global plate motions.

proposed tectonic and dynamic forces acting on the North American plate include gravitational potential energy (GPE) associated with continental surface topography and oceanic ridge push gravity sliding (lithospheric thickening effects), basal tractions from the mantle flow induced by nearby subducting slabs, and boundary shear and normal loads from neighboring plates (Humphreys & Coblentz, 2007). The contribution of lithospheric thickening effects to the global plate motions has been shown to be minor compared with those from subducting slabs (e.g., Lithgow-Bertelloni & Richards, 1998). Ghosh, Becker, et al. (2013) found that accounting for GPE effects slightly increases the predicted plate motion magnitudes over North America but has almost no influence on the directions of the North American plate motion. The effects of the boundary loads are debated due to the complexity of the tectonic blocks near the plate boundaries (e.g., McCaffrey, 2005). The strong North American cratonic root (Laurentia) (Figures 1b and 2) deeply penetrates into the upper mantle (Jordan, 1975), allowing the North American plate to more efficiently couple with mantle flow than the plates with a thinner lithosphere (Conrad & Lithgow-Bertelloni, 2006). It has been shown that basal drag generated by convective mantle flow is able to shift the lower portion of the North American cratonic root westwards along the direction of the plate motion (Kaban et al., 2015). Nevertheless, the primary buoyancy sources contributing to the plate-mantle coupling beneath the North American plate are yet to be investigated. . For each tomography model, only the areas where the seismic velocity is slower than the ambient mantle by more than one standard deviation of the global negative seismic velocity anomalies are counted. The areas with the most votes are confidently considered the slowest seismic velocity zones, including the region around the southwestern North America to the north of the East Pacific Rise (16-18 votes). See Shephard et al. (2017); Hosseini et al. (2018) and the references therein for more detail about vote maps (https://www.earth.ox.ac.uk/∼smachine/cgi/index.php?page=vote_maps) and the seismic tomography models. (b) Seismic velocity perturbations beneath the North American plate at a depth of 150 km from S40RTS tomography models. The continental cratonic root is featured with the fastest seismic velocity, in contrast with the slowest seismic velocity around the southwestern North America within the black dashed frame, which indicates the existence of sublithospheric partial melt.  Table 1), which consists of stiff cratonic keels, modestly stiff slabs, and temperature-dependent viscosity. Viscosity η is the absolute viscosity with unit Pa⋅s. Symbols A, B, and C in the right panel denote the locations of three near-field Aleutian, central American, and Caribbean active subducting slabs around the North American plate, respectively. Other active subducting slab segments are regarded as far-field slabs from the North American plate.
Effects of the enhanced basal tractions from strong plate-mantle coupling to the North American plate motion can be studied via global mantle flow models with lateral viscosity variations (LVVs) incorporating realistic lithospheric plate structures and strong cratonic roots driven by buoyancy estimated from seismic tomography models. A robust mantle flow model should not only predict plate motions, a constraint on horizontal flow; it should also be consistent with geoid or gravity field, which is an important vertical observable constraint on mantle buoyancy and rheology (Forte, 2007;Forte et al., 2015;Ghosh et al., 2010). Only a small number of studies (e.g., Forte, Quéré, et al., 2010;Ghosh, Becker, et al., 2013;Ricard & Bai, 1991) tried to fit both observations in a single dynamically-consistent mantle circulation model. Previous numerical experiments focusing on the North American plate constrained by both plate motions and geoid lead to three major findings. First, the North American plate is coupled to the underlying mantle flow induced by the ancient descending Kula-Farallon plate system in the lower mantle (e.g., Forte et al., 2007;, which results in convergent upper mantle flow beneath the eastern part of the North American continent (e.g., Ghosh, Becker, et al., 2013). Second, including an observationally-constrained, well-resolved global slab geometry model (Gudmundsson & Sambridge, 1998) into upper mantle significantly improves the predicted plate motion over North America (Ghosh, Becker et al., 2013). Ghosh, Becker, et al. (2013) speculate that this improvement in the plate motion is due to the addition of a well-resolved slab beneath the Aleutian Arc. Third, strong sublithospheric divergent flow, driven by the predicted active upwelling below the Pacific margin just off the coast of the southwestern North America, appears underneath the western part of the North American continent Ghosh, Becker et al., 2013). This strong upwelling, potentially rooted in the transition zone and lower mantle , is considerably enhanced by the extremely slow seismic velocity anomalies that are interpreted to be hot and buoyant around southwestern North America in the uppermost mantle ( Figure 1). The strong northeasterly divergent flow at the base of the lithosphere exerts large resisting basal tractions underneath the western U.S. making it hard to fit the North American plate motion (e.g., Ghosh, Becker, et al., 2013).
The first and the second findings attribute different buoyancy sources to the primary driving forces for the North American plate. Is the North American plate mainly driven by the active subducting slabs in the upper mantle (Ghosh, Becker, et al., 2013), or the ancient Kula-Farallon slab remnants in the lower mantle , or both? Moreover, if active subducting slabs in the upper mantle are important, is the North American plate mainly driven by the near-field Aleutian, central American, and Caribbean slabs along the margins (Ghosh, Becker, et al., 2013), or by the large-scale mantle flow induced by other far-field subducting slabs? The contributions to the North American plate motion from the global subduction system require further investigations. The wide distribution of the partial melt around mid-ocean ridges, rifts, and hotspots complicates the interpretation of slow seismic anomalies in the uppermost mantle (Hammond & Humphreys, 2000;Rooney et al., 2012). Using a velocity-to-density scaling (δ ln ρ/δ ln V S ) that assumes seismic velocity can be entirely attributed to temperature will over-estimate the positive buoyancy in regions where partial melt is present. The southwestern North America around the Pacific margin is featured with one of the slowest seismic velocity anomalies beneath the lithosphere (Figure 1), which is a northward extension of the East Pacific Rise under North America (Menard, 1960). The extremely slow seismic wave speed in this region may reflect partial melting in the upper mantle. Previous mantle flow models have considered the lateral perturbations of the velocity-to-density scaling (δ ln ρ/δ ln V S ) arising from cratonic roots (Figure 1b) (e.g., Ghosh, Becker, et al., 2013). However, the potential lateral perturbations of the velocity-to-density scaling associated with the shallow partial melt around the southwestern North America have not been explicitly incorporated into the construction of the density field, which may be over-estimating the positive buoyancy and hence the upwelling beneath this region. The excessively strong divergent flow induced by over-estimated upwelling from the partial melt within uppermost mantle around the southwestern North America may overly resist the North American plate motion, challenging the reliability of the third finding. Thus accounting for the effect of partial melt on seismic velocity could be a significant factor to understand the driving forces for the North American plate.
This study investigates the driving mechanism for the North American plate in a global 3D fully dynamically-consistent instantaneous mantle flow model constrained by both geoid and plate motions. We address the issue on whether the North American plate is mainly driven by near-field marginal subduction zones or far-field subducting slabs, and reexamine the influence of the lower mantle buoyancy sources (including Kula-Farallon slab remnants and the active upwelling below the Pacific margin of the North American plate) on the surface plate motions. We also explore the effect of the lateral perturbations of the velocity-to-density scaling associated with the shallow partial melt around the southwestern North America to the predicted North American plate motion.

Mantle Flow Modeling
We use an open source finite element mantle convection code, ASPECT (Advanced Solver for Problems in Earth ConvecTion) to compute mantle flow and geoid (Bangerth et al., 2019a;Heister et al., 2017;Kronbichler et al., 2012;Liu & King, 2019), which have been well benchmarked in Liu and King (2019). The global instantaneous mantle circulation is modeled in a 3D spherical shell, assuming an incompressible Stokes fluid with Boussinesq approximation with the governing equations for the conservation of mass and momentum: where ⃗ is the velocity, p is the pressure, η is the viscosity, ρ 0 is the reference density, δρ is the density perturbation, ⃗ is the gravity, and , is the strain rate. The reference density ρ 0 and the radial gravity are taken from the PREM profiles (Dziewonski & Anderson, 1981).
The density anomalies are constructed from the S40RTS shear-velocity seismic tomography model (Ritsema et al., 2011), which has been demonstrated to match best with surface deviatoric stress and geoid in a recent comparison study between mantle flow models derived from various seismic tomography models (Wang et al., 2015). The high velocity anomalies underneath cratons as determined from 3SMAC model (Nataf & Ricard, 1996) are removed down to 300 km to make the cratonic roots neutrally buoyant (e.g., Forte & Perry 2000;Jordan, 1978). A constant velocity-to-density scaling of 0.25 is used to derive the density perturbations from S40RTS throughout the mantle (e.g., Ghosh, Becker et al., 2013). We retain the shallow buoyancy from seismic tomography models except for the cratonic roots to include the contribution of lithospheric thickening effects on plate motions. A lateral perturbation to the velocity-to-density scaling associated with the shallow slow seismic anomalies beneath the southwestern North America may be included as discussed later. Additionally, to improve the resolution of the negative buoyancy distribution at narrow subduction zones, a global model of regionalized upper mantle (RUM) slabs inferred from seismicity (Gudmundsson & Sambridge, 1998) is incorporated into the density field in some cases as shown later. To construct this combined buoyancy structure, the density perturbation from RUM slabs is fixed at +40 kg/m 3 , and we eliminate all other fast seismic anomalies in S40RTS from 100 to 660-km depth to avoid an extra contribution from subduction, following the approach in Ghosh, Becker, et al. (2013). We also test four other negative buoyancy values for the RUM slabs at +30 kg/m 3 , +35 kg/m 3 , +45 kg/m 3 , and +50 kg/m 3 , and find that they yield nearly identical results as +40 kg/m 3 .
For the rheology structure, we use a five-layer radial viscosity profile including lithosphere (0-100 km), asthenosphere (100-300 km), upper mantle (300-410 km), transition zone (410-670 km), and lower mantle (670-2891 km), similar to the division in Hager and Richards (1989). Weak plate boundaries defined by the NUVEL-1A plate motion model (DeMets et al., 1994) are incorporated into lithosphere. Strong cratonic keels, extending to 250-km depth, are introduced using the 3SMAC model (Nataf & Ricard, 1996). With the transition zone viscosity fixed at 1 × 10 21 Pa⋅s, we use the combined density field from S40RTS seismic tomography model and RUM slabs to search for the rheology structure that best fits the global geoid and plate motions simultaneously. The best fitting viscosity model has a stiff lithosphere (1.0 × 10 24 Pa⋅s) with weak plate boundaries (1.0 × 10 21 Pa⋅s), weak asthenosphere (8.0 × 10 19 Pa⋅s), moderately weak upper mantle (3.0 × 10 20 Pa⋅s), and moderately strong lower mantle (5.0 × 10 22 Pa⋅s). The viscosity of the most rigid cratonic keels are fixed at 1.0 × 10 26 Pa⋅s, 100 times stronger than the lithosphere (Figure 2a). The temperature-dependent viscosity is incorporated below 100 km depth following (Ratcliff et al., 1996) where η 0 is taken from the five-layer radial viscosity profile, = − ln , is the temperature perturbation from the S40RTS density anomaly (α is the coefficient of thermal expansion), E determines the strength of the temperature-dependence and is set to 0.01, which produces LVVs on two orders of magnitude below lithosphere. Due to the strong positive buoyancy associated with the extremely slow seismic velocity below the southwestern North America, the asthenospheric viscosity below this region from the cases using a constant velocity-to-density scaling of 0.25 is reduced to as low as 10 17 -10 18 Pa⋅s (Figure 2b). These low viscosity magnitudes are consistent with previous local studies using viscoelastic modeling (e.g., Bills et al., 2007). Incorporating the strong cratonic keels, the variation range of LVVs in the asthenosphere extends over eight orders of magnitude ( Figure 2b). We do not test various strength values of temperature-dependent viscosity (E) because (a) the effect of temperature-dependent rheology may be small compared to the variability due to uncertainties in the seismic tomography models ; and (b) aside from modestly increasing the predicted plate velocity, temperature-dependent rheology has almost no effect on the long-wavelength geoid and the orientations of plate motions (Ghosh, Becker, et al., 2013;Moucha et al., 2007).
The strength of subducting slabs remains an open question (Billen, 2008) given that both mantle flow models with strong slabs (e.g., Alisic et al., 2012;Stadler et al., 2010) and weak slabs (e.g., Forte et al., 2015;Ghosh & Holt, 2012) can fit the plate motions well. The fraction of the upper mantle slabs' excess weight that effectively contributes to the slab pull force remains uncertain because various decoupling mechanisms including bending resistance (Chapple & Forsyth, 1979), slab detachment (Spakman et al., 1988), fault zone resistance (van Dinther et al., 2010), and viscous support of slabs (Alpert et al., 2010;Isacks & Molnar, 1971), are hard to accurately constrain. In addition, strong slabs in the mantle flow models using Newtonian rheology without the consideration of plastic yielding in the bending region (Hines & Billen, 2012) or a low-viscosity mantle wedge associated with either non-linear rheology (McAdoo, 1982) or hydration (Billen & Gurnis, 2001) are found to produce geoid lows over subduction zones, contrary to the observations (Moresi & Gurnis, 1996;Yoshida & Nakakuki, 2008). Further considering that the North American plate, majorly as a non-subducting overriding plate, has only a small fraction of marginal subduction zones, it should be primarily driven by the slab suction through basal drag. Therefore, we only include modestly stiff RUM slabs, which are stronger than the ambient mantle by a factor of 3, from 100 to 410 km depth. The slab viscosity is the same as the surrounding mantle from 410 to 660 km depth (e.g., Liu & Stegman, 2011) accounting for the potential weakening effect associated with the grain size reduction as a result of phase transformations in the transition zone (Mohiuddin et al., 2020;Riedel & Karato, 1997). It should be emphasized that our treatment on the slab viscosity cannot rule out the strong slabs in the upper mantle. The issue on the strength of the upper mantle slabs is simply beyond the scope of our Newtonian global circulation models focusing on the North American plate, constrained by both plate motions and geoid.
Our viscosity structure ( Figure 2) is similar to Ghosh, Becker, et al. (2013) except for the lithosphere and cratonic keels. We use a stiffer lithosphere (1.0 × 10 24 Pa⋅s) with weaker plate boundaries (1.0 × 10 21 Pa⋅s) ( Figure 2a). This allows our models to further fit the global plateness constraint (Kreemer et al., 2003) and fully decouple the motions between the surface plates by three orders of magnitude viscosity reduction at weak zones (Kaban et al., 2014). The cratonic keels in our models are uniformly strong (1.0 × 10 26 Pa⋅s) down to 250-km depth (Figure 2b), with no weakening below 100 km in the asthenosphere (see Figure 1 in Ghosh, Becker, et al., 2013). Our stronger deep cratonic keels result in more efficient coupling of basal tractions from mantle flow underneath continental plates, especially the North American plate with one of the largest cratonic roots (Figures 1b and 2b). The 250-km thick strong cratonic keels are underlain by a globally-extending low-viscosity asthenosphere down to 300 km (Dalton et al., 2009;Gung et al., 2003;van Summeren et al., 2012;Yuan & Romanowicz, 2010), given that the deepest less viscous part of the cratonic roots may be sheared and deformed (e.g., Kaban et al., 2015).
In our computations, quadratic finite elements for the velocity and linear finite elements for the pressure (Q 2 Q 1 ) are used (Bangerth et al., 2019b), which satisfies the LBB (Ladyzhenskay-Babuska-Brezzi) stabilization conditions (Brezzi & Douglas, 1988). The whole 3D spherical shell domain is divided into 64 uniform layers with 49,152 cells per layer, leading to the vertical resolution of ∼45 km and the horizontal resolution of ∼0.6 × 0.6°. Because ASPECT uses Q 2 Q 1 finite elements, as compared with Q 1 P 0 (linear velocity and constant pressure) finite elements commonly used in other convection codes, the effective vertical and horizontal resolution could be considered to be doubled (i.e., ∼22.5 km and ∼0.3 × 0.3° respectively). Free-slip boundary conditions are applied to both surface and core-mantle boundary. The pure rotation components for the whole mantle in the flow velocity solutions are removed following Zhong et al. (2008). ASPECT uses consistent boundary flux (CBF) method to determine the dynamic topography contribution to the geoid (Zhong et al., 1993). The predicted geoid is computed in the spectral domain up to spherical harmonic degree 40, including the effects of self-gravitation, and is summed over to the spatial domain (see Liu & King, 2019 for details).

Comparison With Observational Constraints
We use the EGM 2008 geoid model (Pavlis et al., 2008) with the hydrostatic flattening removed following Chambat et al. (2010) as the vertical observational constraint ( Figure 3a). Both cross-correlation (e.g., Becker & Boschi, 2002) and variance reduction (e.g., Ricard et al., 1993) between the model geoid and EGM 2008 geoid are computed from spherical harmonic degrees 2-20 as the contribution from degrees 21-40 is negligible. Cross correlations range from −1 (perfect negative correlation with the observation) to 1 (perfect positive correlation with the observation) and variance reductions are numbers of 1 or less with 1 indicating no reduction in variance relative to the observation. Further details of the two evaluation methods for geoid modeling can be found in Liu and Zhong (2016). We use the NUVEL-1A plate motion model (DeMets et al., 1994) in the no-net-rotation (NNR) reference frame as the horizontal observational constraint ( Figure 3b) to evaluate our predicted plate motions. In accordance with the NNR reference frame, the net rotation of the surface flow from our models is further removed before the comparison considering that lithosphere still has a considerable net rotation over the underlying mantle in the presence of shallow strong LVVs (Becker, 2006;Zhong, 2001). The predicted plate motions are evaluated using three criteria in the spatial domain, direction fit to NNR-NUVEL-1A (ξ), magnitude fit to NNR-NUVEL-1A (β), and plateness (P). The direction and magnitude fits are evaluated by the weighted average (by the product of the two velocity magnitudes) dot product and the average logarithmic magnitude ratio between two velocity vectors respectively (Becker, 2006); the global plateness (0.87) from the NNR-GSRM plate motion model (Kreemer et al., 2003) are used to evaluate the plate-like level of our predicted surface velocity fields considering that there may be certain degree of intraplate deformation of the plate motions derived from seismic tomography-based mantle flow models (Becker, 2006): where ‖ ⋅ ‖ denotes the Euclidean norm and 〈⋅〉 denotes the average calculation. ⃗ is the predicted plate motion vector and ⃗ is the observed NNR-NUVEL-1A plate motion vector. ⃗ = ⃗ |⃗ | and ⃗ = ⃗ |⃗ | , are the normalized unit vectors. ξ (−1 to 1) and β (0: perfect fit; negative: slower prediction; positive: faster prediction) measure the plate direction and magnitude fits, respectively. ⃗ is the velocity vector from the best fitting Euler Pole of the plate polygon S i in the model. is the normalized plate area. P i and P are individual and global plateness, respectively. To ensure the dynamic self-consistency of our mantle flow models, we evaluate our predicted plate motions both globally and over 8 large NUVEL-1A plates (Africa, Australia, Eurasia, South America, North America, Pacific, Nazca, India) individually.  (DeMets et al., 1994) in no-net-rotation (NNR) frame with legend at the top right. The red lines outline NUVEL-1A plate boundaries.

Case 1: Buoyancy From Seismic Tomography Alone
The global geoid and plate motions predicted from considering only S40RTS buoyancy (Case 1 in Table 1) using the best fitting rheology structure described in Section 2.1 are shown in Figures 4a and 4b respectively.
The predicted global geoid is in good agreement with EGM 2008 geoid, with a cross correlation of 0.796 and a variance reduction of 0.623. The peak values (max: ∼110 m; min: ∼−116 m) of the predicted global geoid are also close to those of the observed geoid (Figure 3a). The geoid highs over Africa, central Pacific, and subduction zones, as well as geoid lows around India ocean, Siberia, Hudson Bay, and Antarctic are well reproduced. Overall the pattern of our predicted geoid is comparable to both previous geoid modeling results from only radially-varying viscosity (e.g., Liu & Zhong, 2016;Steinberger & Calderwood, 2006) and the ones considering LVVs (e.g., Ghosh et al., 2010;Yang & Gurnis, 2016). The directions of the plate motions predicted from Case 1 globally fit NNR-NUVEL-1A reasonably well (ξ = 0.858), but the predicted global velocity magnitudes are smaller than the observed magnitudes (β = −0.110). The predicted surface flow is also plate-like with a global plateness of 0.860, which is close to the global plateness of 0.87 from the NNR-GSRM model. However, when using the buoyancy from only S40RTS seismic tomography, the predicted continental plate motions do not fit NNR-NUVEL-1A as well as the oceanic plates do (Table 1; Figures 4b vs. 3b). This is especially true for the North American plate (Table 1), which moves in almost the opposite direction to the NNR-NUVEL-1A model (Figure 5b vs. 5a), resulting in a strong, negative direction fit (ξ = −0.641). The geoid and plate motion results from the S40RTS only buoyancy are consistent with the results from the composite SMEAN seismic tomography model (Becker & Boschi, 2002) using the similar viscosity structure (Ghosh, Becker, et al., 2013). The seismic tomography derived buoyancy alone is again shown to be insufficient to predict continental plate motions, especially the motion of the North American plate. It has been argued that this is because (a) traditional global seismic tomography models typically lack the regional resolution to accurately capture the upper mantle slabs along narrow subduction zones (Lu & Grand, 2016); and (b) the lateral variations in the velocity-to-density scaling are not robustly incorporated into the construction of 3-D mantle buoyancy field to explicitly accommodate thermal and compositional origins respectively, as manifested by the improved fits to surface observable constraints over continental plates (e.g., Forte, Quéré, et al., 2010) via joint seismic-geodynamic inversions for 3-D mantle structures ).

Case 2: Buoyancy From Seismic Tomography Plus Near-Field Slabs
We next add the Aleutian, central American, and Caribbean slab segments from the RUM slab model into the S40RTS buoyancy to better represent the three near-field slabs ( Figure 2b) along the North American plate margins (Case 2 in Table 1). This case, to our knowledge, has not been computed in previous studies, but is pivotal to determine whether the motion of the North American plate can be primarily attributed to the driving forces from near-field subduction zones. With a subtle change to this local negative buoyancy, the model still fits the global geoid with a high cross-correlation (0.809) and variance reduction (0.642) as shown in Figure 4c. The predicted global plate motions (Figure 4d) are also similar to those driven by only S40RTS buoyancy (ξ = 0.872; β = −0.098; P = 0.865). Generally, the predicted continental plate motions still do not fit NNR-NUVEL-1A as well as the oceanic plates do (Table 1; Figure 4d vs. 3b), which is to be expected because we only incorporate the high-resolution negative buoyancy of three local slab segments around the North American plate. Comparing Figures 5c and 5a, the resultant surface velocity directions of the North American plate in Case 2 are still nearly opposite to NNR-NUVEL-1A directions (ξ = −0.612) with only a limited improvement from Case 1 (ξ = −0.641). Clearly, only incorporating the three well-resolved near-field slabs into global buoyancy field is not enough to correct the North American plate motion.

Case 3: Buoyancy From Seismic Tomography Plus Global Slabs
We further insert the remaining far-field slabs from RUM slab model into Case 2 to enhance the resolution of all the upper mantle slab segments (Figure 2b) in the global buoyancy field, which forms Case 3 shown in Table 1. Compared with Case 1 from S40RTS only buoyancy, there is still not much change in the resultant global geoid (Figure 4e), which has a cross-correlation of 0.810 and variance reduction of 0.627. The predicted global plate motions (Figure 4f), on the other hand, are considerably improved with an excellent global direction fit 8 of 21 (ξ = 0.950) and a good global magnitude fit (β = 0.057). With the addition of upper mantle RUM slabs, globally the predicted surface flow, with a plateness of 0.925, is slightly more plate-like than NNR-GSRM model (0.87). The improved global direction fit is mainly attributable to the significant improvement on the orientations of predicted continental plate motions driven by a well-resolved global subduction system in the upper mantle, which also increases the magnitudes of global plate motions (Table 1; Figure 4f vs. 3b). The direction fit of the North American plate motion (ξ = 0.061) as shown in Figure 5d, while still considerably worse than those of the other continental plates (Table 1; Figure 4f), is substantially improved from the poor negative fits in the S40RTS  only buoyancy case (Case 1) and the S40RTS plus near-field RUM slabs case (Case 2). These results are generally consistent with Ghosh, Becker, et al. (2013), who use global RUM slabs in conjunction with the combined seismic tomography model SH_TX2008, which is constructed from a regional S-wave model beneath western U.S. (Schmandt & Humphreys, 2010) nested in the global S-wave model TX2008 ). However, the North American plate motion in our model does not exhibit the rotation pattern around the center of the continent as seen in Ghosh, Becker, et al. (2013). Considering that we use a similar viscosity structure and the same constant velocity-to-density scaling of 0.25 as in Ghosh, Becker, et al. (2013), this subtle inconsistency is presumably due to local differences in the buoyancy structures between S40RTS and combined SH_TX2008 seismic tomography models beneath North America.
The significant discrepancy between the direction fits of the North American plate and the other continental plates in Case 3 indicates that it is not sufficient to incorporate well-resolved global slabs into the buoyancy field to successfully predict the North American plate motion. There remain factors uniquely important for the dynamics of the North American plate that are yet to be accounted for Ghosh, Becker, et al. (2013) found that assigning a slightly larger viscosity, by a factor of ∼10 than other weak plate margins, to the San Andreas Fault plate boundary segment (Platt et al., 2008) further improved the direction fit of the predicted plate motion over North America. However, there is still considerable difference between their predicted and NNR-NUVEL-1A North American plate motions, especially for the western part of the continent. NNR-NUVEL-1A plate motion turns to a southwesterly orientation over this region, roughly perpendicular to the direction of nearby Pacific plate motion ( Figure 5a); whereas their predicted North American plate drifts northwestwards over this region, approximately parallel to the nearby Pacific plate motion (Figure 14c in Ghosh, Becker, et al., 2013). The northwesterly motion of the North American plate at the western part of the continent in their model is potentially generated by the inherent mechanical coupling between Pacific plate and the North American plate as a result of stronger San Andreas Fault. We identify another factor, that is, over-estimation of the positive buoyancy associated with the extremely slow seismic velocity induced by shallow partial melt beneath the southwestern North America, to explain the difficulty to fit the North American plate motion in Case 3.

Cases 4-8: Local Correction of Positive Buoyancy Associated With Partial Melt
If the extremely slow seismic velocity within uppermost mantle around the southwestern North America is the result of widely distributed partial melt, this will lead to over-estimated positive buoyancy based on the velocity-to-density scaling which simply assumes a thermal origin. This will induce exaggerated upward motion and hence excessively strong divergent flow at the base of the lithosphere in this region. To investigate this problem, we locally reduce the velocity-to-density scaling for the negative velocity anomalies within the black dashed frame in Figure 1b in top 300 km, which prescribes a first-order estimation of the extremely slow seismic velocity around the southwestern North America.
To isolate the effect of the over-estimated positive buoyancy around the southwestern North America, we first remove the negative velocity anomalies within this region (δ ln ρ/δ ln V S = 0) to construct the buoyancy field from S40RTS only model (Case 1) to form Case 4, as shown in Table 1 (Table 1). While the direction fit of the predicted North American plate motion from Case 4 ( Figure 5e) is still negative (ξ = −0.230), it is apparently better than that from Case 1 (Figure 5b, ξ = −0.641). The predicted North American plate motion from Case 4 indicates that without incorporating well-resolved global slabs into the buoyancy field, a model with the correction of the positive buoyancy around the southwestern North America alone is not able to explain the motion of the North American plate. Based on Case 4, we further add three near-field Aleutian, central American, and Caribbean slab segments from the RUM slab model (Figure 2b) into the buoyancy field to construct Case 5, which is used to investigate whether the three near-field slabs are sufficient to drive the North American plate in presence of the correction of the local positive buoyancy around the southwestern North America. Globally, the predicted geoid ( Figure 4i) and plate motions (Figure 4j) from Case 5 fit the observations to the similar levels as Case 4 ( Table 1). The direction fit of the predicted North American plate motion from Case 5 (Figure 5f) is 0.396, which is notably better than −0.230 from Case 4. This indicates that the three near-field slabs contribute moderately to the direction of the North American plate motion in presence of the correction of the local positive buoyancy around the southwestern North America. However, the speed of the North American plate in Case 5 (Figure 5f) is extremely slow (β = −0.568). These results further show that the negative buoyancy from only the three near-field slabs are not vigorous enough to drive the North American plate. The large-scale driving effect from other far-field slabs cannot be ignored. Therefore, our next investigation on the strength of the local positive buoyancy around the southwestern North America will be based on the combined density structures constructed from S40RTS and global RUM slabs.
Due to the uncertainty in the geometry of the partial melt (e.g., Hammond & Humphreys, 2000), it's impossible to define a single appropriate velocity-to-density scaling for the strong negative velocity anomalies in this region. Previous joint inversion of global seismic and mantle convection data sets have obtained a radial velocity-to-density scaling profile featured with values less than 0.05 (a factor of 5 smaller than our global constant velocity-to-density scaling value, 0.25) in the top 300 km for non-cratonic ambient mantle (see Figure 3 in . Therefore, we test three different reduced velocity-to-density scaling values at 0.05 (Case 6), 0.02 (Case 7), and 0 (Case 8). Case 6 and Case 8 can be regarded as end-member tests representing the upper and lower limits of the positive buoyancy associated with the effect of partial melting respectively. Case 7 can be regarded as a moderate test representing a reasonable guess of the reduced velocity-to-density scaling for negative velocity anomalies in this region. The results of these three cases are shown in Table 1. With only a subtle local modulation of the global shallow positive buoyancy, the resultant global geoid (Figures 6a, 6c, and 6e) and plate motions (Figures 6b, 6d, and 6f) are nearly the same as those from Case 3 (Table 1). However, when the velocity-to-density scaling value for shallow negative velocity anomalies around the southwestern North America is reduced, the direction fits of the predicted North American plate motions (Figures 7a, 7b, and 7c) are significantly improved from Case 3 (Table 1). Without the southwestern shallow positive buoyancy (Case 8), the predicted North American plate motion (Figure 7c) matches the NNR-NUVEL-1A model best (ξ = 0.753).
Reduced velocity-to-density scaling values of 0.05 (Case 6) and 0.02 (Case 7) also yield similar good direction fits of the predicted North American plate motions (Figures 7a and 7b, respectively) to the NNR-NUVEL-1A model (Case 6: ξ = 0.734; Case 7: ξ = 0.746). Indeed, considering the lateral perturbation of the velocity-to-density scaling associated with the shallow partial melt around the southwestern North America significantly improves the predicted North American plate motion. We are aware that the reduced velocity-to-density scaling also increases the local asthenospheric viscosity around the southwestern North America based on the temperature-dependent rheology as shown in Equation 3, which may increase the mechanical coupling between the Pacific and the North American plate, and hence improve the direction fit of the predicted North American plate. However, we have verified that the effect of the increased asthenospheric viscosity to the North American plate motion is negligible compared with the significant effect of the reduced positive buoyancy in this region (see Supporting Information S1). We also note that the best predicted North American plate motions in our models (Case 6: β = −0.182; Case 7: β = −0.176; Case 8: β = −0.172) are slower than that in the NNR-NUVEL-1A model. We attribute this to the lack of GPE effect associated with the continental surface topography over North America in our models (Ghosh, Becker, et al., 2013).

Case 9: Exclusion of Lower Mantle Buoyancy
We finally test one more model to remove the S40RTS buoyancy below 660 km depth in Case 7 to evaluate the importance of the lower mantle buoyancy (Case 9 in Table 1). Compared with Case 7, the global geoid ( Figure 6g) fit from Case 9 is significantly degraded with a cross correlation of 0.624 and a variance reduction of 0.379. This is not surprising because the lower mantle buoyancy has been shown to contribute a significant power to the long-wavelength geoid (Hager & Richards, 1989;Richards & Hager, 1988). Without the lower mantle buoyancy, each plate still fits the NNR-NUVEL-1A directions very well, although the motions become systematically slower (Table 1). Correspondingly, the global plate motions (Figure 6h) from Case 9 maintain an excellent direction fit of 0.968 to the NNR-NUVEL-1A model with an overall slightly slower magnitude (β = −0.071). The global plateness (0.920) is also nearly the same as that from Case 7. The predicted North American plate motion (Figure 7d) yields an excellent direction fit of 0.915 (even better than 0.746 in Case 7) to NNR-NUVEL-1A with smaller velocity magnitudes (β = −0.290). Aside from decreasing the speed of the predicted surface velocity, exclusion of the lower mantle buoyancy does not considerably impact the patterns of global plate motions. Ghosh, Becker, et al. (2013) speculates that the near-field slab beneath Aleutian Arc may play an important role to improve the predicted North American plate motion with the incorporation of high-resolution upper mantle subducting slabs (Gudmundsson & Sambridge, 1998). However, the limited improvement on the predicted North American plate motion with the addition of the well-resolved near-field Aleutian, central American, and Caribbean slabs (Case 1 vs. Case 2; Case 4 vs. Case 5) demonstrates that near-field slabs are not the primary negative buoyancy sources that drive the North American plate. Furthermore, the substantially improved fit of the North American plate motion with the addition of well-resolved global upper mantle slabs (Case 2 vs. Case 3; Case 5  Table 1), (c, d) S40RTS with a reduced velocity-to-density scaling of 0.02 for the negative velocity anomalies around the southwestern North America plus global RUM slabs (Case 7 in Table 1), (e, f) S40RTS with the negative velocity anomalies around the southwestern North America removed plus global RUM slabs (Case 8 in Table 1), and (g, h) only the upper mantle (top 660 km) S40RTS with a reduced velocity-todensity scaling of 0.02 for the negative velocity anomalies around the southwestern North America plus global RUM slabs (Case 9 in Table 1), respectively, using the best fitting rheology structure described in Section 2.1. The red lines outline NUVEL-1A plate boundaries.

Far-Field Slabs Dominantly Driving the North American Plate
vs. Case 8) indicates that far-field subducting slabs play an non-negligible role in driving the North American plate. The active subducting slabs around the western Pacific margins constitute the major part of the global subduction system in the upper mantle (Figure 2b), and hence should provide the primary far-field driving forces. The large-scale mantle flow induced by slabs on the other side of the Earth exert basal tractions that act to drag the North American plate westwards. The coupling of these far-field slab suction forces is enhanced by the strong North American cratonic root (Figures 1b and 2b). The sublithospheric flow patterns of Case 1, Case 2, and Case 3 at a depth of 181 km beneath the North American plate are shown in Figures 8a, 8b, and 8c, respectively. The stiff North American cratonic root moves northeasterly in almost the opposite direction to the NNR-NUVEL-1A model using seismic tomography only buoyancy (Case 1). When the well-resolved three near-field slabs are incorporated into the buoyancy field (Case 2), the northeasterly motion of the cratonic keel persists. However, the trend of northwesterly motion of the cratonic keel begins to emerge when other well-resolved far-field slabs are further incorporated into the buoyancy field (Case 3). The large-scale mantle flow induced by far-field slabs exerts the dominant basal tractions to drag the North American plate.
In addition to the North American plate, other continental plates also strongly couple to the large-scale mantle flow driven by far-field slabs through their rigid cratonic keels, as demonstrated by the improved direction fits of the predicted African, Australian, Eurasian, and South American plate motions using a combined buoyancy field from global seismic tomography and RUM slabs (Case 3). The importance of the large-scale suction due to  Table 1), (b) S40RTS with a reduced velocity-to-density scaling of 0.02 for the negative velocity anomalies around the southwestern North America plus global RUM slabs (Case 7 in Table 1), (c) S40RTS with the negative velocity anomalies around the southwestern North America removed plus global RUM slabs (Case 8 in Table 1), and (d) only the upper mantle (top 660 km) S40RTS with a reduced velocity-todensity scaling of 0.02 for the negative velocity anomalies around the southwestern North America plus global RUM slabs (Case 9 in Table 1), respectively, using the best fitting rheology structure described in Section 2.1. The red lines outline NUVEL-1A plate boundaries. The gray area denotes the North American craton.
far-field slabs in the continental plate motions also explains why oceanic plate motions are well predicted using seismic tomography only buoyancy (Case 1), while continental plate motions can only be well predicted using a combined buoyancy field from both seismic tomography and RUM slabs. The large-scale negative buoyancy from far-field slabs does not appear to be well accounted for in seismic tomography only buoyancy due to the limited resolution in subduction zones, as manifested in the relatively poor direction fits of the continental plates in Case 1. The oceanic plates, on the other hand, are primarily driven by the suction and pull forces from the dominantly one-sided, near-field subducting slabs. The accuracy of the slabs imaged by seismic tomography appears to be sufficient to represent these near-field slab suction forces. When well-resolved global RUM slabs are combined with the seismic tomography model (Case 3, Cases 6-8), the far-field slab suction forces are represented well enough to exert considerable shear tractions at the base of cratonic keels, which significantly improves the predicted continental plate motions; whereas the directions of oceanic plate motions remain nearly unchanged with only the increased speed due to the stronger suction forces from higher resolution near-field slabs (Table 1).
Global negatively buoyant slabs are dominant density sources driving both oceanic and continental plates, but their dynamic mechanisms are different in terms of the strong coupling between large-scale mantle flow driven by far-field slabs and deeply penetrating rigid cratonic roots. This far-field slab suction applied to continental plates is physically consistent with the classical theory of the surface response functions from deep mantle buoyancy sources (Hager & O'Connell, 1981). The importance of far-field slabs to drive continental plate motions, especially the North American plate motion, can be regarded as an extension from the previous theory of nearby slab  Table 1), (b) S40RTS plus near-field Aleutian, central American, and Caribbean RUM slabs (Case 2 in Table 1), (c) S40RTS plus global RUM slabs (Case 3 in Table 1), and (d) S40RTS with a reduced velocity-to-density scaling of 0.05 for the negative velocity anomalies around the southwestern North America plus global RUM slabs (Case 6 in Table 1), respectively, using the best fitting rheology structure described in Section 2.1. The purple dots denote the North American cratonic root. A saturation of the colors occurs when the amplitude of the vertical flow exceeds ±5 cm/yr. The sublithospheric flow patterns beneath the North American plate calculated from Case 7 and Case 8 in Table 1 are similar to that from Case 6. suction forces to drive overriding continental plates (e.g., Conrad & Lithgow-Bertelloni, 2002). Plate tectonics is better regarded as a self-organizing interacting system, as opposed to an isolated balance of forces on each plate. Regional geodynamic modeling focusing on continental plates must account for the effects of large-scale mantle flow induced by far-field subduction zones.

Buoyancy From Shallow Slow Seismic Velocity Associated With Partial Melt
The reduction effect of partial melt on the seismic velocity could be much stronger than that on the density, which may result in the over-estimation of the positive buoyancy from the shallow hot mantle including partial melt based on the velocity-to-density scaling assuming purely thermal origin. Global mantle convection models incorporating the process of partial melting constrained by plate motion history indicate recent considerable melt production at mid-ocean ridges around Gulf of California (Li et al., 2016), as manifested by the extremely slow seismic velocity in this region (Figure 1). Previous mantle flow models neglecting the lateral perturbations of the velocity-to-density scaling with the assumption that these slow seismic anomalies are entirely due to temperature overpredict the positive dynamic topography over the southwestern North America (Ghosh, Becker, et al., 2013). This indicates that the positive buoyancy in the uppermost mantle beneath this region is over-estimated because dynamic topography kernels are most sensitive to shallow buoyancy sources. It is also noteworthy that the radial velocity-to-density scaling profiles from a joint inversion of global seismic and mantle convection data sets are featured with significantly small values below 0.05 at top ∼300 km for non-cratonic ambient mantle (see Figure  3 in . This in turn implies the wide distribution of the partial melt featured with extremely slow seismic velocity in the uppermost mantle, which requires considerable reduced velocity-to-density scaling coefficients to avoid over-estimation of the shallow positive buoyancy. Our numerical tests on the significantly reduced velocity-to-density scaling within uppermost mantle (top 300 km) around the southwestern North America (Cases 4-8) confirm that a reasonable estimation on the positive buoyancy in this region featured with widely distributed partial melt is critically important to explain the dynamics of the North American plate. We further illustrate the importance of this local buoyancy correction by comparing the sublithospheric flow beneath this region between Case 3 without a reduced velocity-to-density scaling and Case 6 with a reduced velocity-to-density scaling of 0.05 (Figure 8c vs. 8d). The northeasterly divergent flow patterns around the southwestern North America are found in both of the two models but the magnitude of the latter one is substantially reduced, which decreases the resisting basal tractions underneath the North American plate. Correspondingly, the stiff North American cratonic root in Case 6 ( Figure 8d) exhibits more robust northwesterly motion than that in Case 3 ( Figure 8c). There is a clear convergent sublithospheric flow pattern beneath the eastern part of the North American continent along the edges of the cratonic root in Case 6 (Figure 8d), consistent with previous findings (e.g., Ghosh, Becker, et al., 2013). It is worthwhile to emphasize that the upwelling and the corresponding sublithospheric divergent flow beneath the southwestern North America still exist in our models with a reduced velocity-to-density scaling in the uppermost mantle (Figure 8d), as these signals are largely driven by the deeper-seated buoyant mantle upwelling (depths >400 km) below this region in the transition zone and the lower mantle . However, without the correction of the underlying shallow positive buoyancy from the extremely slow seismic velocity associated with the partial melt, this northeasterly sublithospheric flow and the resulting resisting basal tractions beneath the North American plate are significantly over-estimated (Figures 8a, 8b, and 8c), which makes the predicted North American plate motion hard to reconcile with the observation. It is also notable that our models (Cases 4-8) with a constant reduced velocity-to-density scaling coefficient within uppermost mantle around the southwestern North America are simply based on a first-order estimation on the presence of the shallow partial melt within a rectangular area (Figures 1b and 8d). We conjecture that a better-defined distribution of the positive buoyancy (e.g., depth dependence and lateral variations) around the southwestern North America could further improve the North American plate motion fit. However, seismic velocity perturbations are jointly determined by temperature and melt in this region. Additionally, the effect of melt on seismic velocity reduction complexly depends on both the melt geometries and the quantity of melt (Hammond & Humphreys, 2000). Therefore, rather than arbitrarily searching for an optimal model by trial and error, we believe that a more complete understanding of the distribution of partial melt and the thermal structure within uppermost mantle around the southwestern North America is necessary.
The shallow partial melt is also widely distributed in other locations within uppermost mantle, such as Eastern Africa (Figure 1a), which contains both the young African-Arabian large igneous province (LIP) and the active East African Rift (EAR) system (Rooney et al., 2012). In contrast with the North American plate, the predicted motions of other major plates are not significantly impacted by neglecting the correction of the shallow positive buoyancy from extremely slow seismic velocity associated with the partial melt in other locations (Case 3 in Table 1). This indicates that the influence of the over-estimated buoyancy from the presence of partial melt within uppermost mantle on the predicted plate motions depends on their relative distributions with respect to the nearby plates. However, the detailed geodynamic investigations into other plates still need to address the shallow positive buoyancy within partial melt areas to robustly estimate the local sublithospheric flow and the corresponding regional surface observables, such as dynamic topography and gravity.
The estimation of the positive buoyancy associated with the presence of the partial melt in the uppermost mantle has been overlooked in previous global mantle flow modeling. Steiner and Conrad (2007) performs an overall study on the effect of the upwelling induced by the slow seismic velocity anomalies in global tomography and concludes that the positive buoyancy associated with slow seismic velocity anomalies generally degrades the fit to global plate motions. However, they employ a constant velocity-to-density scaling assuming a thermal origin throughout the mantle and their discussion focuses on the slow velocity material in the lower mantle. Global mantle flow models based on joint seismic-geodynamic inversions for 3-D mantle density structures incorporating thermal and compositional heterogeneity  show that interpreting the slow seismic velocity anomalies in the lower mantle as thermo-chemical origin improves the fit to plate motions and geoid Simmons et al., 2007). Rowley et al. (2016) further show that the strong upwelling driven by the positive buoyancy in the middle and lower mantle below the East Pacific Rise may considerably contribute to plate driving forces for Nazca and Pacific motions.
Although the implications of the slow seismic velocity in the uppermost mantle has not been explicitly investigated previously, we can divide the past dynamically consistent global mantle flow models into four categories based on the fit to the North American plate motion and verify the consistency between our new finding on the over-estimated positive buoyancy associated with the presence of the partial melt and the previous results. The first group of models are based on the slab only negative buoyancy constrained from history of subduction and they predict the correct pattern of the North American plate motion (e.g., Becker & O'Connell, 2001, Conrad & Lithgow-Bertelloni, 2002Lithgow-Bertelloni & Richards, 1998;Zhong, 2001). This is not surprising because the slab only negative buoyancy does not contain any positive buoyancy induced by the negative density anomalies that resists the North American plate motion. The second group of models are constructed by merging the upper mantle subducting slabs and the lower mantle seismic tomography (e.g., Alisic et al., 2012;Ghosh & Holt, 2012;Ghosh, Holt et al., 2013;Stadler et al., 2010;Wen & Anderson, 1997). Without the positive buoyancy in the upper mantle, the predicted North American plate motions from these models also fit the observation very well. The third group of models are based on the global seismic tomography removing the velocity anomalies in the uppermost mantle within top 200 ∼300 km (e.g., Becker & O'Connell, 2001;van Summeren et al., 2012).
Although the main purpose of the exclusion of the shallow velocity perturbations is to account for the effect of the tectosphere where compositional differences are likely to cancel out the fast anomalies beneath cratons (e.g., Forte & Perry, 2000;Jordan, 1978), the extremely slow seismic velocity anomalies in the uppermost mantle are removed together. Therefore, these models still achieve a moderate good fit to the North American plate motion. The forth group of models are based on the global seismic tomography, potentially combined with the well-resolved upper mantle subducting slabs, excluding the shallow (top ∼300 km) fast seismic velocity anomalies beneath the cratons (e.g., Becker, 2006Becker, , 2017Ghosh, Becker, et al., 2013). These models poorly match the North American plate motion due to the over-estimated shallow positive buoyancy in the regions with the distribution of the partial melt, including the southwestern North America.
The first three groups of models are constructed from incomplete mantle buoyancy structures and hence they are actually not Earth-like models. Lacking the positive buoyancy in the upper mantle, these models do not include the contribution of lithospheric thickening effects, may not produce correct asthenospheric flow patterns, and cannot fit the surface observables that are sensitive to the shallow buoyancy signals (e.g., dynamic topography). Our investigations on the strength of the sublithospheric upwelling around the southwestern North America based on the prediction of the North American plate motion indicate that the shallow buoyancy structures in the forth group of global mantle flow models can be retained to build more Earth-like models as long as one can employ an appropriate velocity-to-density scaling in the regions featured with extremely slow seismic velocity associated with the distribution of the partial melt.

The Role of Lower Mantle Buoyancy
It has been proposed that the dominant basal tractions driving the motion of the North American plate are exerted by the convective flow associated with the descent of the ancient Kula-Farallon plates in the lower mantle below the North American continent (Bokelmann, 2002). This local deep downwelling and an active upwelling below the Pacific margin of the North American plate within the lower mantle are further investigated by , who argue that these two buoyancy sources have a dominantly strong signature in the surface plate motion over North America. Our test on the exclusion of the lower mantle (below 660 km) buoyancy from seismic tomography (Case 9) excellently predicts the directions of the plate motion over North America, but the velocity magnitudes are considerably decreased from Case 7 (Table 1; Figure 7d vs. 7b). We do not overinterpret the best direction fit of the predicted North American plate motion when removing the lower mantle buoyancy because better-defined buoyancy and distribution of the partial melt around the southwestern North America may further moderately improve the direction fit of the North American plate motion in Case 7 to a similar level (see Section 4.2). However, the decrease of the velocity magnitudes over North America is robust, which is seen in the results of other plates as well (Case 9 in Table 1; Figure 6h). The local deep downwelling induced by the ancient descending Kula-Farallon plates and the active upwelling below the Pacific margin in the lower mantle primarily influence the magnitudes, instead of the directions, of the North American plate motion. Globally, lower mantle buoyancy only increases the amplitudes of the surface plate motions (Case 9 in Table 1; Figure 6h). The patterns of the global plate motions are predominantly driven by the active subducting slabs within the upper mantle in terms of suction and pull forces, as discussed in Section 4.1.

Conclusions
In this study, we construct dynamically consistent 3D global instantaneous flow models simultaneously constrained by both geoid and plate motions to investigate the driving mechanism for the North American plate, which is featured with a strong cratonic keel and only a small fraction of marginal subduction zones.
We find that the near-field Aleutian, central American, and Caribbean slabs only have a minor suction effect to drive the North American plate, different from the speculation in Ghosh, Becker, et al. (2013); whereas the large-scale mantle flow induced by other far-field slabs, primarily the major segments around western Pacific subduction margins, provides the dominant driving forces for the North American plate through the coupling to its rigid cratonic keel. These far-field driving forces are also important for the motions of other continental plates with strong cratonic roots. Given that remote subduction generates large-scale flow that drives continental plate motions, it is critical to account for these far-field buoyancy forces in the regional geodynamic investigations over continental plates. The proposed far-field slab suction effect, which is an extension of the previous nearfield slab suction theory, suggests a new form of active plate interactions within the global self-organizing plate tectonic system. We further verify that the over-estimation of the shallow positive buoyancy associated with the partial melt around the southwestern North America results in excessively strong northeasterly sublithospheric divergent flow, which overly resists the North American plate motion. To avoid impairing the ability of global mantle flow models to explain the dynamics of the North American plate, a significantly reduced velocity-to-density scaling (0.05 or smaller in our models) for the negative seismic shear-velocity anomalies within uppermost mantle around the southwestern North America must be incorporated into the construction of the buoyancy field. Correctly addressing the positive buoyancy within the areas including widely distributed partial melt in the uppermost mantle is potentially important for the regional tectonic modeling of other plates as well. Future high-resolution mantle circulation modeling constrained by various surface observables could consider incorporating an accurately-defined global partial melt buoyancy and distribution model into the construction of buoyancy field within uppermost mantle to more robustly estimate regional-scale expressions of mantle flow.
Lower mantle buoyancy, which includes the local deep downwelling induced by the ancient descending Kula-Farallon plates and the active upwelling below the Pacific margin of the North American plate, influences the amplitudes of both North American and global plate motions, but does not control the patterns of the surface plate motions.