Secular change of true polar wander over the past billion years

The rate of movement of Earth’s solid shell relative to its spin axis, or true polar wander, depends on variations in mantle convection and viscosity. We report paleomagnetic and geochronologic data from South China that constrain the rate of rapid true polar wander (>5° per million years) between 832 million years and 821 million years ago. Analysis of the paleomagnetic database demonstrates secular change of true polar wander related to mantle cooling and thermal structure across supercontinent cycles. True polar wander rates are relatively muted with a partially insulated mantle during supercontinent assembly and accelerate as mantle thermal mixing reestablishes with supercontinent breakup. Decreasing true polar wander rate through the Neoproterozoic was succeeded by overall smaller variations in the Phanerozoic. We propose that Neoproterozoic extensive plate tectonic activities enhanced mantle cooling, giving rise to a reduction in mantle convective forcing, an increase in mantle viscosity, and a decrease in true polar wander rates into the Phanerozoic.

Additional description of geochronologic results A medium-grained gabbro-diabase sample (TS01C), composed primarily of plagioclase and pyroxene, was selected for dating. 206 Pb/ 238 U vs. 207 Pb/ 235 U ratios for the analyzed five zircons of TS01C are plotted on the concordia curve ( fig. S3A). Filled ellipses and vertical bars are zircon grains included in the age calculation ( fig. S3). The weighted mean 206 Pb/ 238 U age calculated by the four zircon analyses is 831.51 ± 0.32 Ma, giving the crystallization age of the diabase sample.
Rock magnetic results Thermal magnetic susceptibility experiments for representative samples reveal a significant decrease in the susceptibility around 580°C, suggesting magnetite is the major magnetic carrier of the samples from the sills and country rocks ( fig. S4A-G). An increase in susceptibility after ~425°C in samples and the cooling curves not resembling the heating curves indicate magnetic mineralogy changes after heating experiments. Thermal demagnetization data show no abnormal change in direction or magnitude of the remanence in the high-temperature range ( fig. S5), showing that the potential alterations of magnetic grains do not influence the remanent magnetization. Fractional magnetization vs. temperature diagram ( fig. S4H) manifests two substantial intensity losses between 500-580°C and 300-350°C. The former corresponds to the unblocking temperature range of fine-grained (titano) magnetite; the latter overlaps with the typical unblocking temperature of pyrrhotite at ~320°C and is only observed for samples from the HGS and JZP sections. The 500-580°C magnetization drop appears to be less remarkable for samples from the JZP section, a result of relatively lower intensity held in the high temperature range compared to the total magnetization ( fig. S4H). Together, these results suggest that magnetite is the main magnetic carrier of the remanence. Pyrrhotite is inferred to be present in part of the samples that accounts for the sudden intensity decline at ~320°C.   Table S2. (H) Fractional magnetization vs. thermal demagnetization for representative samples. Between two demagnetization steps, the intensity loss is calculated as the magnitude of the demagnetized vector, instead of the arithmetic difference between the remanent magnetizations. Baked-baked country rocks of the sill JZP01. Unbaked-unbaked country rocks in the JZP section. Vertical grey bands highlight the temperature interval of 300-350°C and 500-580°C. Directions are shown in geographic coordinates. Each title signifies the name of the rock unit and the specimen (Table S2). Vectors in colors indicate the least-squares fits for each remanence component (60). NRM-natural remanent magnetization; LT1 and LT2-low-temperature components; HT1-high-temperature component carried by sills and baked rocks; HT2-hightemperature component carried by unbaked rocks in the JZP section.
Low temperature components (LT1 and LT2) Most of the analyzed samples, including sills and country rocks, carry LT1. LT1 was typically removed from room temperature to ~450°C during thermal demagnetization, unless intermitted by the removal of MT between ~300-350°C for part of the samples ( fig. S5). The mean of LT calculated from all samples in geographic coordinates is D = 357.5°, I = 47.3°, n = 193, α95 = 1.9°, close to the local present field direction and fails a fold test (fig. S6). These suggest LT1 is a viscous remanent magnetization (VRM) of the local present geomagnetic field.
LT2 resides in samples from the Huguosi (HGS) and Jingzhanping (JZP) sections but not from the Heiwan (HW) section. LT2 was isolated from a relatively narrow range, between ~300-350°C, entailing the typical unblocking temperature of pyrrhotite (~320°C) ( fig. S5). LT2 has variable magnetization that accounts for about 10-40% of the samples' total natural remanent magnetization (NRM) (fig. S4H). The mean of LT2 of all samples in geographic coordinates is D = 171.8°, I = -43.9°, n = 135, k = 47.7, α95 = 1.8° ( fig. S6). In the vector component diagrams where LT2 is present, the removal of LT2 between ~300-350°C displays a visible directional shift, showing an opposite trend to the LT1 removal after which LT1 resumes decay. LT2 also fails a fold test and is interpreted as secondary ( fig. S6). However, the age and origin are unclear. Both unstable components, LT1 and LT2, are not further evaluated in the study.
High-temperature component (HT1) and positive regional fold tests We average the specimen directions of each sill to calculate the site-mean directions of HT1 of the seven sills. In geographic coordinates, HT1 display scattered distributions that fall into three groups-a NW-up direction of sills from the HGS and JZP sections, a NW-down direction of HW01, and a NE-up direction of HW02 ( fig. S7). After tilt-correction using the beddings of each sill's country rocks (sedimentary rocks of the lower Fanjingshan Gp), HT1 became significantly clustered (D = 351.3°, I = -44.1°, N = 7, k = 46.5, α95 = 8.9°) ( fig. S7). The seven site-mean directions pass the McFadden (23) fold test at the 95% confidence level [critical ξ = 3.1; ξ1 (insitu) = 3.8; ξ1 (tilt-corrected) = 1.4]. A progressive unfolding test (24) indicates that the highest precision parameter (k) is obtained at 98.9% unfolding ( fig. S7), enveloped by 89-110% suggested for an optimal positive fold test (73). Together, the positive regional fold tests suggest HT1 was acquired before the tilting and folding of the Fanjingshan Gp when the sills were parallel to country rocks of the lower Fanjingshan Gp.
We also compute the HT1 mean using all the specimen directions (specimen mean) to test the stability of HT1 directions calculated by different statistical approaches. The specimen mean of  S7). Compared with the specimen mean, the mean of the site-mean directions is a better representative of HT1 because, in this way, it avoids the potential bias toward the groups of specimen directions with larger sample sizes.  (Table S2). Color-filled/colorrimmed squares represent lower/upper hemisphere mean directions; stars indicate the means; ellipses show radius of 95% confidence cone of the means (α95). Stars indicate the means of the seven site-mean directions. Dec-declination; Inc-inclination; k-precision parameter; Nnumber of directions to calculate the mean; HGS-Huguosi section; HW-Heiwan section; JZP-Jinzhanping section. (C) Progressive unfolding test (24) for HT1 using the site-mean directions of the seven sills. The red vertical line at 98.9% unfolding indicates where k reaches its maximum. The dashed vertical lines bound 90-110% unfolding preferred for a positive fold test suggested by (73). Positive baked-contact test on HT1 In the JZP section ( fig. S1, fig. S9), sites TR53-54-73, baked siltstones collected from 1-15 cm away from the contact of JZP01, yield a stable high-temperature component (480-550°C) similar to HT1 of JZP01 ( fig. S9, fig. S13), interpreted to be a partial thermal remanent magnetization (pTRM) imprinted by the sill. Part of the baked samples shows non-origin-trending decay during ~480-550°C (without change of direction), implying the possibility of a weak residual magnetization associated with the siltstones > 550°C. Nevertheless, fitting great circles to these specimens (between ~480-550°C) does not find another consistent direction near the unblocking temperature (~580°C) but the planes all commonly intersect at HT1 ( fig. S10), indicating such a residue, if present in part of the specimens >550°C, is scattered and unstable. Therefore, the line fits between ~480-550°C do not include the origin for pure isolation of the thermal imprint by the sill. We performed a Watson's V statistic test (74), and the null hypothesis of a common mean for the high-temperature components from JZP01 and the baked siltstones cannot be rejected at 95% confidence level (Watson V = 2.5, smaller than the critical value of 6.9).
Samples from sites TR55-56, unbaked siltstones about 30-200 cm below the intrusive contact of the JZP01, do not preserve a stable and consistent high-temperature component. Sites TR87-88, unbaked grey-green siltstones between JZP01-JZP04, yield a stable HT2 that is significantly different from HT1 ( fig. S9, fig. S13I). HT2, isolated from ~500-580°C and likely held by magnetite, is a unique NE-up direction that does not resemble any other component found from the Fanjingshan region. It is unclear whether HT2 is a primary detrital remanent magnetization or not without a paleomagnetic confidence test. Collectively, the above results demonstrate a positive baked contact test for HT1, indicating it is a primary thermal remanent magnetization acquired during the cooling of the sills.

The Fanjingshan pole and paleosecular variation test
The seven site-mean directions of the sills in sedimentary coordinates are first converted to corresponding VGPs; averaging the seven VGPs gives the Fanjingshan pole (34.7°S, 118.2°E, A95 = 8.6) (Fig. 1A, fig. S11A). To test if our paleomagnetic results have sufficiently recorded paleosecular variations of geomagnetic fields, we analyze the VGP angular standard deviation, SB, after (75) and use the n-1 Jackknife resampling method to obtain the confidence intervals of the mean (76). After eliminating random errors introduced by the sampling and measuring process, SB, calculated for the total 7 VGPs, has a mean of 10 (21,64,67,68). This implies the mafic magmatism in this area might have already ceased by the upper Fanjingshan Gp deposition when the bedding of the lower Fanjingshan Gp remained horizontal. (ii) No unconformity or depositional hiatus is found within the Fanjingshan Gp and the adjacent formations show conformable contact (21,64,67), suggesting no regional tectonic movements that have tilted the strata of the Fanjingshan Gp prior to regional deformation associated with the Jiangnan Orogen. Together with the positive fold test on HT1, this precludes a possibility that the sills emplaced during the country rocks were tilted rather than horizontal. (iii) The U-Pb LA-ICP-MS 822 ± 8 Ma and ca. 835 Ma youngest populations of detrital zircon (64,77) and the 832.0 ± 8.5 Ma tuff from the upper Fanjingshan Gp (69) indicate the depositional age of the upper Fanjingshan Gp may be as young as <832 Ma and imply the emplacement of the ca. 832 Ma sills were likely contemporaneous with part of the Fanjingshan Gp deposition. The evidence above supports the inferred temporal sequence that the mafic sills emplaced before the regional deformation of the Fanjingshan Gp, validating the tilt-correction method and the directions of HT1.     Table S3 and S5. Coloration according to pole ages is applied to the Neoproterozoic poles; Phanerozoic poles are shown in grey.
Summary of TPW records since the late Mesoproterozoic Secular TPW is now processing on Earth at ~1.1° Myr -1 primarily due to contributions of ongoing glacial isostatic adjustment, quantified by the geodetic observations of relative speed of platehotspot motions (79). Recent numerical predictions on the present-day TPW rate, which adopt a more accurate ice age rotational theory, suggest a composite nature of the observed rate resulted from both the residual effects of ice age and modern global sea-level rise (80).
In the last ~100 Myr, Earth's TPW has been considered to operate in a confined manner of only <6°. The relatively muted TPW magnitude and rate can be ascribed to high mantle viscosity (81), a stable and flattened Earth figure introduced by a steady pattern of plate subduction and largescale mantle upwelling, or remanent bulge stabilization/Earth's excess ellipticity (1,81,82). Recently, paleomagnetic data from Scaglia Rossa limestone of Italy provides evidence for a ~12° TPW oscillation from 86 to 78 Ma, with a mean rate of ~3.0° Myr -1 (9). The rapid TPW rate was attributed to a larger convective forcing than the normal (9). This interpretation is consistent with our predictions for invigorated mantle convective driving that accelerated TPW following the mantle thermal remixing with Pangaea breakup after ca. 170 Ma. The amplitude of TPW was estimated to be 12 ± 3° for single trip, corresponding to 24 ± 6° for the cumulative TPW roundtrip. Thereby, we estimate a mean rate of 3 ± 0.75° Myr -1 for the 86-78 Ma TPW oscillations.
Constraints on TPW is more ambiguous in the late Paleozoic. From 400 Ma to 250 Ma, ~60° coherent northward migration of Pangaea at a steady rate of 0.4° Myr -1 from the southern latitudes has been recently suggested as evidence for TPW rather than lithosphere-mantle motion (39). Considering a ~10° uncertainty in the paleolatitude for Pangaea at 400 Ma and 250 Ma [ Fig. 2 from (39)], we estimate the rate uncertainty of ± 0.1° Myr -1 . The authors also found that the mean TPW rate could be larger (0.5-0.6° Myr -1 ) using their updated Pangaea reconstruction (39). Combining the two sources of uncertainties, we estimate the credible intervals of the 400-250 Ma TPW rate to be 0.3-0.7° Myr -1 , with the mean rate at 0.5° Myr -1 .
However, uncertainties remain in the inferred magnitude and rate of TPW due to the inability to distinguish tectonic-motion components from the total apparent polar wander path of Gondwana, during the closure of the Rheic Ocean and amalgamation between Gondwana and Laurussia (7). Accordingly, Torsvik et al. (7) considered a maximum TPW rate at ~2.0° Myr -1 at 430-410 Ma interpreted for a synthesized apparent polar wander path for Laurentia, concerning a portion of the pole shift was allocated to plate-tectonic movement. Taking their reported great-circle difference (GCD) as an approximation to the pole uncertainty and using the standard uncertainty propagation, we quantify the means and uncertainties of the maximum TPW rates using their original data [ Table. 5 from Torsvik et al. (7)] that yield 1.9 ± 0.7° Myr -1 between 430 and 420 Ma and 2.0 ± 0.6° Myr -1 between 420 and 410 Ma.
Late Ediacaran to Cambrian TPW records (555-505 Ma) have been proposed from Laurentia and Gondwana (87), among which the most recent robust data set came from Australia (88). According to updated paleomagnetic results from the Amadeus Basin, central Australia, the paleogeographic shift of assembled Gondwana was ~60° over ~40 Myr, with plate velocity estimates from 8-28 cm yr -1 , inferred to be either nonuniformitarian plate tectonics or an episode of rapid TPW, or both (88). The Early Cambrian TPW hypothesis has been interpreted as consistent with coeval global d 13 Ccarb variations and exceptionally fast evolutionary rates in the biosphere (a factor of 20 higher), both of which were potentially triggered by prominent paleogeographic shifts of global continents and sea-level fluctuations due to large-amplitude and long-term TPW (87,89).
Earlier in the Ediacaran, ca. 615-565 Ma rapid and, potentially, oscillatory TPW has been postulated based on paleomagnetic data from both sedimentary and igneous rocks in Laurentia, West Gondwana, Baltica, and Australia (14). A ~90° CW rotation of these continents between ca. 615 and 590 Ma is followed by an opposite ~90° rotation around a similar TPW axis from 575 to 565 Ma (14). From global data, TPW mean rate estimates for the two TPW phases were recently proposed at ~3 ± 1° Myr -1 and ~10° Myr -1 , respectively (14). Notably, the ~10° Myr -1 estimate for the latter 575-565 Ma event comes from averaging very scattered results from various continents, which also have very large uncertainties (error bars from +/-2° Myr -1 to +/-15° Myr -1 ) (14). Robert et al. (15) evaluated Ediacaran TPW mechanisms with global geographic reconstructions and a mantle dynamics model. They found that the reactivation of a girdle of subduction surrounding the continents followed by a reduced return flow could potentially explain the observed TPW oscillation. The modeled mean TPW velocities are consistent between the two TPW shifts: ~2.9° Myr -1 and 3.2-3.5° Myr -1 for 615-590 Ma and 575-565 Ma, respectively. For the 575-565 Ma event, although the mean TPW rate estimate from observations (~10° Myr -1 ) is larger than the modeled mean values (3.2-3.5° Myr -1 ), they nevertheless overlap each other within uncertainty. Therefore, we regard the modeled speed of 3.2-3.5° Myr -1 from (15) as a reasonable and likely estimate, while ~10° Myr -1 is a possibility of a maximum speed of this interval. Alternatively, these large Ediacaran pole shifts have been interpreted to record non-Geocentric Axial Dipole (GAD) fields with extremely weak geomagnetic moments and rapid reversals (90,91).
In the Tonian, during the Bitter Springs isotopic stage (810-795 Ma), three consecutive paleomagnetic poles from the carbonate units of the Akademikerbreen Group, East Svalbard, Norway, exhibit two nearly opposite shifts in the direction of ~83° total rotation within ~15 Myr (10). These shifts in pole position broadly coincide with changes in d 13 Ccarb and local sea level, supporting a TPW interpretation (10). The age of the three Svalbard paleomagnetic poles that bracket the Bitter Springs isotopic stage were constrained by chemostratigraphic correlations with the sedimentary sessions in Canada and Ethiopia, in which absolute dating were obtained (11,92,93). Using these age constraints, the Svalbard poles yield a TPW rate of ≥40-54 cm yr -1 (≥3.6-4.9° Myr -1 ) for the whole Bitter Springs Stage (from ca. 810 to <795 Ma) (11). Using a Monte Carlo approach accounting for the uncertainties between pole ages, Park et al. (12) determine the TPW rate to be ≥60 cm yr -1 (≥5.4° Myr -1 ) between the pre-Bitter Springs pole (> ca. 810 Ma) and syn-Bitter Springs pole (ca. 810-795 Ma). From South China, the newly reported Xiajiang pole (816-810 Ma, pre-Bitter Springs isotopic stage) (12) and Madiyi pole (ca. 805 Ma, syn-Bitter Springs isotopic stage) (32) constrain pole motion significantly less than that predicted from the coeval Svalbard poles, pointing to the possibility of TPW rotation counteracted by tectonic motion of South China.
Preceding the Bitter Springs Stage, the ca. 832 Ma Fanjingshan pole from this study and the ca. 821 Ma Xiaofeng pole (28,29) constrain the total pole motion of South China to be 54.7 ± 13.9° at an estimated rate of 3.5-7.8° Myr -1 resulting from pure TPW (accounting for possible tectonic movements). If using the proposed revised Xiaofeng pole from Jing et al. (94), the total pole motion between 832 and 821 Ma would be 69.8 ± 16.9° at a rate of 7.1 ± 1.7° Myr -1 , agreeing with the results using the original Xiaofeng pole (28) within uncertainty. This large low-to-high-latitude motion of South China at 832-821 Ma is comparable to the ~56° pole-to-equator latitude change of Baltica in an overlapping period (Fig. 1C-D), determined by the 848 ± 27 Ma ( 40 Ar-39 Ar biotite plateau date) (35)/834 ± 9 Ma (Rb-Sr isochron age) (95) Hunnedalen dykes pole (35) and ca. 836-803 Ma Katav Formation pole and Inzer Formation pole (16,17) loosely constrained by Rb-Sr illite isochron dates from shale in the Inzer Formation (96). As Baltica has been established as part of the assembled Rodinia in the Tonian (10,36), the coherent motion of Baltica (assembled Rodinia) and possibly disconnected South China reinforce the hypothesis that these motions recorded TPW.
Before the Neoproterozoic, it has been suggested that TPW may account for, at least in part, the rapid plate velocity (~27-34 cm yr -1 ) implied by the Keweenawan Track between ca. 1,110 and 1,080 Ma obtained from the Midcontinent Rift, Laurentia (56). With the implementation of the Bayesian approach, Swanson-Hysell et al. (56) have evaluated different models to disentangle plate-tectonic and TPW components of the 1,110-1,080 Ma apparent polar wander path. The TPW rate could be feasibly interpreted to be 5-22 cm yr -1 accompanied by a fast plate-tectonic motion of Laurentia at ~15 cm yr -1 . Alternatively, the Keweenawan Track can potentially be explained entirely as plate tectonic motion, with TPW close to zero (56). Similar rapid pole motions of this stage are supported by possibly contemporaneous paleomagnetic results from carbonate successions of the Nanfen Formation, North China, that show significant pole-toequator latitudinal change consistent with the data from Laurentia (97). Nevertheless, the timing of these pole shifts is poorly constrained without direct and precise age controls on these carbonate units (97).   Table S8. For the analysis, observations are randomly selected from the original data set (with replacement) to form 1,000 resampled data sets of the same sample size. For each re-sampled data set, the best-fit slope is determined and plotted (grey line). The three cases (A), (C), and (E) demonstrate a set of sensitivity tests: (A) assuming non-GAD origins of the rapid pole shifts during the Ediacaran and therefore no estimates of TPW of this time; (B) assuming the likely and conservative Ediacaran TPW rate estimates consistent with those modeled by Robert et al. (15); (C) assuming the possibility of maximum TPW speeds (~10° Myr -1 between 575 and 565 Ma) implied from Robert et al. (14). TPW rate (° Myr -1 ) is translated to the corresponding plate velocity (cm yr -1 ) at a point 90° from the TPW axis using a factor of 11.1. . For all the three cases, the 95% confidence intervals of the slope estimates overlap each other and do not envelop slope = 0. These results suggest that the positive correlation of TPW rate vs. age (long-term decrease of TPW rate with time) is (i) significant at 95% confidence level and (ii) independent from the variably interpreted rapid TPW rates or their presence in the Ediacaran. It shows that most cases (~67%) place the simulated rates above 10 cm yr -1 , weighted toward the endmember of very high tectonic velocities. Therefore, the distribution adequately accommodates large effects of potential tectonic motions. Notes : Colored fractions indicate zircons used in age calculation. X-internal (analytical) uncertainty in the absence of all external or systematic errors; Y-incorporates the U-Pb tracer calibration error; Z-includes X and Y, as well as the uranium decay constant errors. MSWD-mean square of weighted deviates; n -number of analyses included in age calculation; N -total number of zircon analyses.    Table S3.

Neoproterozoic paleomagnetic poles for South China
Notes: The 0.6 notation-paleomagnetic poles after f = 0.6 inclination correction; ! and "-latitude and longitude of paleomagnetic poles in present-day geographic coordinates; dp/dm -semiminor and semimajor axes of 95% polar error ellipse; A 95 -radius of 95% confidence cone of paleomagnetic pole.   Notes: ! and "-latitude and longitude of paleomagnetic pole in present-day geographic coordinates; dp/dm-semiminor and  Notes : (1) Considering no plate tectonic movements. (2) Accounting for plate tectonic movements. Data file S1. Fits of high-temperature components. Dg/Ig: declination and inclination in geographic coordinates. MAD: maximum angle of deviation. Fitting method signifies fits forced/unforced to pass through the origin. Fitting steps correspond to the demagnetizing steps as recorded in the measurement-level data files.