Embracing Uncertainty to Resolve Polar Wander: A Case Study of Cenozoic North America

Our understanding of Earth's paleogeography relies heavily on paleomagnetic apparent polar wander paths (APWPs), which represent the time‐dependent position of Earth's spin axis relative to a given block of lithosphere. However, conventional approaches to APWP construction have significant limitations. First, the paleomagnetic record contains substantial noise that is not integrated into APWPs. Second, parametric assumptions are adopted to represent spatial and temporal uncertainties even where the underlying data do not conform to the assumed distributions. The consequences of these limitations remain largely unknown. Here, we address these challenges with a bottom‐up Monte Carlo uncertainty propagation scheme that operates on site‐level paleomagnetic data. To demonstrate our methodology, we present an extensive compilation of site‐level Cenozoic paleomagnetic data from North America, which we use to generate a high‐resolution APWP. Our results demonstrate that even in the presence of substantial noise, polar wandering can be assessed with unprecedented temporal and spatial resolution.

resolve the history of that plate's latitudinal and rotational motion, for as long as the record remains sufficiently dense.
Paleomagnetic data are typically assembled and reported in a hierarchical fashion (Tauxe et al., 2010). First, individual paleomagnetic sample directions are grouped and averaged into sites. Sites are geological units for which the magnetic remanence direction of all associated samples is interpreted to have recorded the same spot reading of the geomagnetic field. Given that the timescale of remanence acquisition is typically quite short relative to the timescale of secular variation of Earth's magnetic field, individual sites do not average out this variability. A virtual geomagnetic pole (VGP) calculated from site-level data thus represents the position of the magnetic dipole axis at the time of magnetization acquisition. Traditionally, the subsequent step is to average site-level VGPs to yield a study-level "paleomagnetic pole" (hereafter "paleopole"), which, in the framework of the time-averaged geocentric axial dipole hypothesis, is interpreted to correspond to Earth's spin axis. Finally, study-level paleopoles are conventionally aggregated to construct an APWP representing the apparent motion of the spin-axis (e.g., Besse & Courtillot, 2002;Kent & Irving, 2010;Torsvik et al., 2012). Note, that the additional step of first averaging to paleopoles in this framework, rather than averaging VGPs from all locations, has no logical basis (Hansma & Tohver, 2020), but the simplicity of working with paleopoles is attractive and has led to wide adoption. Commonly, APWPs are generated by means of moving average techniques, wherein a Fisher mean (R. Fisher, 1953) is computed from all the paleopoles whose age falls within a running window (Irving, 1977). An alternative approach employs a spherical spline (Jupp & Kent, 1987) to fit smooth curves to assemblages of paleopoles (e.g., Torsvik et al., 2008).
These traditional approaches to synthesize paleomagnetic data into APWPs have evolved little in the last decades, but they include several deficiencies: they do not incorporate spatial and temporal uncertainty, and assume all paleopoles to contain equal weight and to conform to a Fisher distribution function (R. Fisher, 1953). Despite increasing attention to the fact that uncertainties are not routinely propagated through the paleomagnetic hierarchy (e.g., Heslop & Roberts, 2016), the incorporation of site-level uncertainties into APWPs remains problematic because it requires a non-linear transformation from direction-to pole-space that renders traditional parametric assumptions (i.e., Fisher statistics) on the pole level ineffective (Heslop & Roberts, 2020). Recently, Vaes et al. (2022) showed that constructing APWPs from VGP-level data instead of from paleopoles enables knowledge of uncertainties to be incorporated in the hierarchy. However, paleomagnetic data are commonly compiled at the study-level and generating compilations of site-level data can be an onerous task, so they took the approach of simulating pseudo-VGPs based on the Fisher statistics of paleopoles. The propagation of uncertainties from deeper hierarchical levels (e.g., site-level and measurement-level data) remains an open problem.
Here, we introduce a new approach for synthesizing APWPs that focuses on data compiled at the site-level. We overcome historical limitations using actual data and its associated uncertainty to implement a bottom-up Monte Carlo uncertainty propagation scheme that avoids making assumptions about distribution statistics on the pole level.

Data and Methods
Constructing APWPs from site-level data requires novel tools and data compilations. Historically, paleomagnetic data have been reported at the site-level but compiled at the study-level. Here, we compiled Cenozoic (0-60 Ma) site-level data and associated metadata from North America (Section 2.1). In Section 2.2, we present a new framework enabling APWPs to be constructed from such data.

Cenozoic North America at the Site-Level: Compilation and Structure
Site-level paleomagnetic data include a wealth of information that is lost when collapsed into summary statistics describing a paleopole. Our compilation retains this information. Here we note a few salient aspects of the compilation, while the full details are presented in Text S1 in Supporting Information S1. First, we retrieve individual site-level paleomagnetic directions and their corresponding VGPs. This approach enables us to work with the actual observational assemblages of directions/VGPs, rather than parametric models of them. Second, each site-level datum is associated with a numeric age and uncertainty that may not be represented in the overall study-level paleopole. Each site-level datum is thus liberated from its study-level collection, enabling it to be independently re-combined with coeval site-level data. Third, we preserve a comprehensive set of metadata associated with each site-level datum that enables filtering by a variety of arbitrary thresholds (e.g., number of samples, within-site scatter). Finally, with site-level data, state-of-the-art paleomagnetic statistical tests (e.g.,  can be applied to the data sets associated with each pole. We focused our site-level case study on Cenozoic North America (0-60 Ma), given the large volume of data previously published. Our compilation, which includes 1,796 site-level records extracted from 39 studies (Figure 1), is available within the MagIC database (http://dx.doi.org/10.7288/V4/MAGIC/19673). Further descriptions of each study are provided in Text S1 in Supporting Information S1.

Monte Carlo Uncertainty Propagation Scheme
We consider the problem of fitting an APWP from a temporal sequence of instantaneous geomagnetic field readings (i.e., site-mean directions) from the same tectonic plate. The implied position of the magnetic pole at the time of magnetization acquisition (VGP) can be calculated from each site-mean direction, and its observed distribution reflects two key components: (a) scatter arising from secular variation of the magnetic field and (b) plate motion relative to the spin axis. Given that secular variation at all latitudes produces VGPs that are circularly symmetric (Deenen et al., 2011;Tauxe & Kent, 2004), the moving Fisher average technique can be used to estimate the "time-averaged" field for each age, that is assumed to be a geocentric axial dipole.
We need an approach to propagate the spatial and temporal uncertainty associated with each input datum to the estimated APWP. To incorporate these uncertainties, estimates were generated by Monte Carlo resampling of the site-mean directions and their associated ages coupled with weighted bootstrap resampling of these site-level data ( Figure 2). This provides a means to propagate uncertainties from the site-level into the APWP by simultaneously resampling both the data set and directions. Our workflow has only mild assumptions: directional errors are assumed to be Fisher-distributed around the site-mean directions (e.g., Lanos et al., 2005) prior to their transformation into VGPs. Additionally, site-level ages are considered to have normally or uniformly distributed uncertainties, depending on whether they are directly radiometrically dated or stratigraphically constrained, respectively (e.g., Rose et al., 2022;Thébault & Gallet, 2010).
When VGPs are distributed unevenly in a window, a mean pole position calculated through a uniformly weighted moving window could be biased toward positions associated with a different age than the mean window ageparticularly when there is a high temporal density of observations to one side of the window. To address this bias, we utilized a weighted bootstrap resampling approach, where we decrease the probability of including VGPs within each window in inverse proportion to their distance from the central age. Moreover, we assign the age of the window as the mean age of the distribution of the resampled VGPs. This approach reduces artifacts associated with uneven VGP age distributions, making the mean position more representative of the center of the window.
We summarize our algorithm in the following steps (see Data Availability Statement for implementation). We use all the site-level data from the collection of study-level paleopoles, then, 1. We generate a simulated sample of the complete site-level data set via Monte Carlo resampling. This involves generating resamples of the site-mean direction in the original data set using a Fisher distribution with a mean equal to the original site-mean direction and dispersion parameter. We take the same number of samples as in the original site mean and calculate the Fisher mean to determine the resampled site-mean direction. If the original age is constrained by a single radiometric date, we sample a new age from a normal distribution with mean equal to the original age and reported uncertainty. If it is stratigraphically constrained, we sample from a uniform distribution with given lower and upper bounds. 2. Each Monte Carlo resampled site-mean direction is then transformed to a VGP. The result is a simulated VGP position with an assigned age for each site-level datum. 3. We then compute the moving Fisher mean with a given window-width parameter. However, the uneven distribution of VGP ages in a window increases the bias of the mean pole position toward ages associated with a higher temporal density of observations than the mean window age. This effect is mitigated by sampling each VGP from a non-uniform probability distribution during the bootstrap step. We assign a sampling probability proportional to 1 − |t i − t 0 |/2(t max − t min ) for each VGP, where t i is the age of the VGP, and t 0 , t max , and t min represent the central, maximum, and minimum ages of the moving window, respectively. After sampling the VGPs from each window, we compute the Fisher mean pole position and assign an age value based on the mean of the age distribution of the VGPs (rounded to the nearest integer). This method is equivalent to running a Fisher running mean with a weighted moving window. 4. We then calculate the implied rate of apparent polar wander (APW) and the shape of the distribution of VGPs, defined from the eigenvalues ratio of the orientation matrix (Scheidegger, 1965). Additionally, we record the number of independent VGPs involved in each moving window, as well as the number of distinct studies from which these VGPs are obtained.
By repeating this procedure (e.g., 10 4 times), we generate empirical distribution functions of statistically plausible mean pole positions for each age. A nonparametric estimate of the time-varying mean pole is computed as the principal eigenvector of the orientation matrix of the statistical ensemble (e.g., Scheidegger, 1965;Tauxe et al., 1991). Furthermore, the confidence region can be evaluated through the spread of vectors around the mean pole (Gallo et al., 2018). The quantile 95% of the empirical distribution of great circle distances from the estimated mean pole provides a circular confidence region (Θ 95 , Heslop & Roberts, 2016). In addition, by excluding 5% of the most distant vectors, confidence limits can be constructed independently for latitude and longitude. Alternatively, a kernel-density estimate can be used to draw a nonparametric 95% confidence region (e.g., N. Fisher & Hall, 1989;Gallo et al., 2017).
It is important to note that the shape of the APWP is dependent on the size of the moving window, which should be large enough to average out secular variation without obscuring polar wander. In other words, small window widths result in an overfitted model (i.e., a curve that contains artifacts due to under-averaged secular variation) while large windows result in underfitting of the model, obscuring polar wander due to plate motion. We find the optimal balance between overfitting and underfitting (i.e., window width optimization) of the moving Fisher mean by defining a stability metric that quantifies variations in the APWP between different splits of the data set (Text S2 in Supporting Information S1). The optimal value of the window width can then be found using the elbow technique of Satopaa et al. (2011).

A High-Resolution APWP for Cenozoic North America
Our new approach propagates uncertainties from the site-level into the final APWP, and selects the optimal value for the width of the moving (Fisher) average in a fully data-driven way. Using this approach, we present a 1-Myr resolution APWP for Cenozoic North America (Figure 3, Figure S2 and Table S1 in Supporting Information S1), available in a Github repository associated with this paper (https://github.com/PolarWandering/Gallo_ etal_2023_APWP_construction/blob/main/data/NAM_Cenozoic_APWP.csv). Although our database consists of 1,796 records, we based our analysis on the original authors' filtering of the sites, that is, 952 VGPs derived from 39 studies. Future work considering alternative quality selection criteria could broaden our understanding of their effects (e.g., Gerritsen et al., 2022).
Figures 3b-3d decomposes the APWP into pole latitude, pole longitude, and APW rate time-series to enable a more straightforward assessment of the path. These time series all smoothly vary over the past 60 Ma, but after 28 Ma there is a distinct and progressive increase in North America's paleolatitude throughout the remainder of the time-series. This corresponds with a linear segment of our new APWP, following a tighter loop in the APWP from ca. 60 to 28 Ma. Interestingly, despite the changes to the APWP trajectory at ca. 28 Ma, the rate of APW remains rather constant (within uncertainty) throughout the considered time interval.
We compare APWPs resulting from moving averages on paleopoles (e.g., Torsvik et al., 2012), moving averages on simulated VGPs , and moving averages with a weighted window on the actual VGPs (this study) in Figure 4a. In Figure 4d, we illustrate the angular distance between them. The angular distance between the classic moving average path and the pseudo-VGP generated path is smaller than the distance between the former and the path from our new approach. This behavior is expected given that the pseudo-VGPs are sampled from the probability distribution functions of the paleopoles used in the classic moving average and both paths are generated with a uniform moving window. Overall, all three paths have similar characteristic features. However, the new path shows a more progressive rate of change from 28 Ma to the present-day, in contrast to the increased APW rate (more pronounced latitudinal change) between ca. 20 and 10 Ma in the other paths. This difference results from the application of a uniform running mean to an uneven distribution of poles/VGPs through this time interval by the prior approaches. Without compensatory weighting, the lower data density between 10 and 20 Ma (Figure 1c) skews the pole positions toward data at the edges of the moving windows, leading to an artifact of an apparently faster rate.

Advancements in APWP Construction
Traditional approaches of APWP construction have been instrumental in delineating first-order polar wander histories. However, these approaches have limitations that may affect APWP accuracy by: (a) introducing an artificial hierarchical level, (b) not propagating uncertainties across hierarchical levels, and (c) assuming probability distribution functions that may be inadequate to model the data. Crucially, how these shortcomings manifest as bias in APWPs, or distort their reported uncertainties, remains largely unexplored. Vaes et al. (2022) provided insight into the consequences associated with (a) by reporting a correlation between the number of VGPs in a paleopole and its precision-as measured by its distance from the reference pole to which it contributed. Thus, by ascribing each paleopole equal weight, a mean pole could be biased toward paleopoles with insufficient information to average out secular variation. Although N-dependent weighting schemes can be implemented to seek to compensate for this bias (e.g., McFadden & McElhinny, 1995;Wu et al., 2021), a more robust solution is to directly compute polar wander from VGPs, preventing the introduction of artificial hierarchies. The construction of APWPs directly from assemblages of instantaneous readings of the paleomagnetic field (e.g., VGPs) moreover aligns with the fundamental hypothesis that the time-averaged field represents a geocentric axial dipole (McElhinny et al., 1996). Recent studies have taken such an approach by constructing APWPs directly from VGPs (Hansma & Tohver, 2020) or by sampling the assumed probability distribution functions of paleopoles (Vaes et al., , 2023. The construction of APWPs from VGPs (or pseudo-VGPs) also innately propagates study-level uncertainties into an APWP (Hansma & Tohver, 2020;Vaes et al., 2022Vaes et al., , 2023 thereby making progress toward addressing limitation (b). Our framework takes this approach further, enabling uncertainties to be propagated from the site-level.
With modern data-sharing initiatives, our algorithm can be straightforwardly extended to propagate sample-or measurement-level uncertainties (Heslop & Roberts, 2016;Khokhlov & Hulot, 2016).  Information S1). In these instances, the classic approach to APWP construction presents two options: dismiss the data, or impose Fisher statistics regardless of the underlying data. In both cases, information is lost, and in the latter case propagation assuming a Fisher distribution (such as in Vaes et al. (2022)) may introduce further bias. In our framework, there is no need to impose probability distribution functions on paleopole or VGP distributions. Following theoretical expectations of a time-averaged field, we use Fisher statistics to estimate a measure of centrality, and through repeated calculations we build empirical distribution functions whose shapes can be monitored. This flexibility also provides a more robust way to incorporate uncertainties that otherwise require the adoption of alternative parametric models (e.g., Pierce et al., 2022).
Discussions of uncertainty in polar wander have often overlooked the importance of the temporal dimension, with conventional APWP approaches treating the (estimated) magnetization ages of paleopoles as an unambiguous age without incorporating uncertainty. The failure to incorporate temporal uncertainties can lead to significant bias. Prior efforts toward temporal uncertainty propagation (e.g., Gallo et al., 2021;Hansma & Tohver, 2020;Swanson-Hysell et al., 2014;Vaes et al., 2023) represent advances but worked from the study-level which has important limitations. First, the relative temporal relationships between sites in a given study are collapsed. Second, applying a single parametric distribution at the study level may be a poor representation of the distribution of site ages and does not enable incorporation of additional geochronology that may be available at the site level. Our framework overcomes these problems and associated biases by enabling age uncertainties to be defined at, and propagated from, the site-level, in the form of both numerical age assignments and the relative stratigraphic ordering of sites.
Our approach provides additional important advantages. Input data filters and quality selection criteria (e.g., number of samples, within-site scatter) heavily influence the outcome of paleomagnetic studies (Gerritsen et al., 2022), and implicit in all APWP construction approaches is the assumption that the base paleomagnetic data have been processed in the same way. In practice, this is almost never true. For many compilations, the retroactive homogenization of the base data requires enormous effort, and any sensitivity testing is prohibitively time-consuming. Because our framework operates on compilations with a comprehensive set of meta-data, it allows the standardized treatment of the base data (from the site-level) via the application of any arbitrary selection criteria.

Geodynamic Applications
The ability to generate high-resolution APWPs presents exciting new opportunities for detailed cross-comparisons between plate kinematic reference frames. Although there is no direct functional relationship between APW and absolute plate motion rates because paleomagnetic data do not directly constrain longitude, APWPs manifest lithospheric motions (Gordon et al., 1979) associated with both changes in latitude and plate orientation that can provide geodynamic insights. Differences between paleomagnetic and hotspot-based reference frames can reveal rotations of the entire solid Earth relative to the spin-axis-a phenomenon referred to as true polar wander (Goldreich & Toomre, 1969). A richer comparison between these reference frames would necessitate the recovery of full plate motion (i.e., including longitude) from paleomagnetic data, which has long been considered an alluring but ultimately intractable problem. However, recent research (e.g., Gallo et al., 2022;Rose et al., 2022) has reinvigorated paleomagnetic Euler pole analysis (Gordon et al., 1984), and shown that it can retrieve paleo-kinematic information from assemblies of paleomagnetic data of sufficiently high-resolution. While conventionally constructed APWPs cannot meet these requirements, our new framework may provide APWPs that can. Future applications of paleomagnetic Euler pole to this new high-resolution APWP for North America are warranted. Other outstanding geodynamics questions elsewhere could be tackled through the generation of additional high-resolution APWPs for other continents, and in deeper time.

Outlook
Two significant obstacles to the broad application of our methodology remain: first, the compilation of site-level data beyond that presented here will require considerable effort, requiring community coordination. The increasing utilization of shared data archives (e.g., MagIC, Jarboe et al., 2012) encourages community coordination associated with the compilation of paleomagnetic data at all hierarchical levels. A second, more daunting challenge is presented by the heterogeneity and general sparsity of paleomagnetic data across space and time. This study focuses on Cenozoic data from North America because it presents one of the most data-rich pockets in space and time. Most other spatiotemporal windows have far fewer data (although our compilation also reveals a noteworthy dearth of paleomagnetic data from stable North America between 10 and 20 Ma). In particular, we would echo Vaes et al. (2023) in calling attention to the need for additional paleomagnetic data from stable continental interiors, which, despite prevailing opinion to the contrary, are still far too scarce.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
All data are included in a Github repository associated with this paper (https://github.com/PolarWandering/Gallo_ etal_2023_APWP_construction) that is also archived in Zenodo (https://doi.org/10.5281/zenodo.7860336). This repository provides an accessible implementation of this analysis for community use via Jupyter notebooks (Kluyver et al., 2016). Additional analysis was performed using PmagPy (Tauxe et al., 2016), Scipy , and NumPy (Harris et al., 2020). Graphics were created using seaborn (Waskom, 2021) Heslop). We are grateful to Trond Torsvik, Carmen Gaina, and the Centre for Earth Evolution and Dynamics (CEED) for supporting this early career project via the YoungCEED programme. We also thank Lei Wu for discussions in support of this work. We appreciate thoughtful reviews of the manuscript from Lisa Tauxe and Andy Biggin as well as the editorial handling of Monica Korte. Special thanks to Aurora, whose timely arrival added fresh perspectives and enriched the development of this work.