Bayesian Paleomagnetic Euler Pole Inversion for Paleogeographic Reconstruction and Analysis

Apparent polar wander paths (APWPs) synthesized from paleomagnetic poles provide the most direct data for reconstructing past paleogeography and plate motions for times earlier than ca. 200 Ma. In this contribution, we describe a new method for APWP synthesis that extends the paleomagnetic Euler pole analysis of Gordon et al. (1984, https://doi.org/10.1029/TC003i005p00499) by placing it within the framework of a Bayesian inverse problem. This approach incorporates uncertainties in pole positions and age that are often ignored in standard treatments. The paleomagnetic Euler poles resulting from the inversions provide estimates for full‐vector plate motion (both latitude and longitude) and associated uncertainty. The method allows for inverting for one or more Euler poles with the timing of changepoints being solved as part of the inversion. In addition, the method allows the incorporation of true polar wander rotations, thus providing an avenue for probabilistic partitioning of plate tectonic motion and true polar wander based on paleomagnetic poles. We show example inversions on synthetic data to demonstrate the method's capabilities. We illustrate application of the method to Cenozoic Australia paleomagnetic poles which can be compared to independent plate reconstructions. A two‐Euler pole inversion for the Australian record recovers northward acceleration of Australia in the Eocene with rates that are consistent with plate reconstructions. We also apply the method to constrain rapid rates of motion for cratonic North America associated with the Keweenawan Track of late Mesoproterozoic paleomagnetic poles. The application of Markov chain Monte Carlo methods to estimate paleomagnetic Euler poles can open new directions in quantitative paleogeography.


Introduction
Plate tectonics is the differential motion of near-rigid blocks of lithosphere across the surface of Earth, separated by relatively narrow regions of deformation in spreading centers, transform faults, and subduction zones. The relative rigidity of plates means that the motion of most of Earth's surface can be described by a set of Euler poles which specify the position (in latitude and longitude) of a rotation axis and a rate of rotation about this axis for a given plate (cf. Cox & Hart, 2009). Individual points on a plate undergoing rigid rotation are described by small circle paths (Figure 1).
Euler poles are widely used for describing plate motions due to their simplicity and compactness (e.g., Argus et al., 2011;DeMets et al., 2010). Motion between tectonic plates has a tendency to remain constant, or approximately so, over time scales of millions of years (e.g., Iaffaldano et al., 2012;Müller et al., 2016). This consistency of motion can be seen physically expressed in the shape of oceanic fracture zones and in hotspot tracks across the lithosphere. These features often form gently curving arcs over large portions of Earth's surface that can be well Abstract Apparent polar wander paths (APWPs) synthesized from paleomagnetic poles provide the most direct data for reconstructing past paleogeography and plate motions for times earlier than ca. 200 Ma. In this contribution, we describe a new method for APWP synthesis that extends the paleomagnetic Euler pole analysis of Gordon et al. (1984, https://doi.org/10.1029/TC003i005p00499) by placing it within the framework of a Bayesian inverse problem. This approach incorporates uncertainties in pole positions and age that are often ignored in standard treatments. The paleomagnetic Euler poles resulting from the inversions provide estimates for full-vector plate motion (both latitude and longitude) and associated uncertainty. The method allows for inverting for one or more Euler poles with the timing of changepoints being solved as part of the inversion. In addition, the method allows the incorporation of true polar wander rotations, thus providing an avenue for probabilistic partitioning of plate tectonic motion and true polar wander based on paleomagnetic poles. We show example inversions on synthetic data to demonstrate the method's capabilities. We illustrate application of the method to Cenozoic Australia paleomagnetic poles which can be compared to independent plate reconstructions. A two-Euler pole inversion for the Australian record recovers northward acceleration of Australia in the Eocene with rates that are consistent with plate reconstructions. We also apply the method to constrain rapid rates of motion for cratonic North America associated with the Keweenawan Track of late Mesoproterozoic paleomagnetic poles. The application of Markov chain Monte Carlo methods to estimate paleomagnetic Euler poles can open new directions in quantitative paleogeography.
Plain Language Summary Movement of Earth's tectonic plates results in large changes in continent positions over time known as paleogeography. Reconstructing paleogeography is central for understanding long term changes to Earth's surface and interior. Ocean floor data are very useful for determining these motions, but only exist for the most recent 4% of Earth's history. For older times, we use magnetic directions in rocks that show positions relative to the north pole. These records are called paleomagnetic poles and plate motions result in progressions of paleomagnetic poles called paths. Earth scientists need ways to synthesize poles into paths. Ideally, such a synthesis uses uncertainties on pole positions and ages. We developed a new way to reconstruct plate motion from paleomagnetic poles using an approach called Bayesian inversion. We apply the method to pretend data to illustrate its capabilities. We then apply the method to real data such as that from Australia over the past 60 million years. This analysis shows change in the continent's position and the rate of its motion that is similar to approaches that rely on ocean floor data. For older times, the method can help reconstruct positions and provide estimates of plate motion rates and associated uncertainty.

ROSE ET AL.
The stage pole description of plate motions has been widely used in both continental reconstructions (e.g., Boyden et al., 2011) and in geodynamical modeling (e.g., Bull et al., 2014;McNamara & Zhong, 2005). Most reconstructions of plate motions over the past 200 million years rely heavily on fitting Euler pole rotations to oceanic fracture zones, hotspot tracks, seafloor magnetic isochrons, and, to a lesser extent, paleomagnetic data (Muller et al., 1993;Müller et al., 2016;Seton et al., 2012). These relative plate motions are used to construct plate circuits that relate plates to one another and to the spin axis (e.g., Müller et al., 2016;Torsvik & Cocks, 2017). However, as we look further back in Earth history, many of the records on which these plate tectonic reconstructions rely largely disappear due to the subduction of oceanic lithosphere. Given the lack of oceanic crust older than ∼200 Ma, the paleomagnetic record from continental rocks is the dominant remaining data. These paleomagnetic data can be used in conjunction with geological data that provides information on the tectonic setting of continental margins as well as additional information, such as the correlation of geologic terranes, to develop paleogeographic reconstructions.
It is more challenging to reconstruct past plate motions from paleomagnetic pole positions than from data derived from oceanic lithosphere for a number of reasons, including: (a) the data are often sparser; (b) traditional paleomagnetic analysis constrains paleolatitude and the orientation of a continental block, but does not constrain paleolongitude without additional assumptions; and (c) paleomagnetic poles can have appreciable uncertainty in both position and age. Gordon et al. (1984) noted that progressions of paleomagnetic poles (which are commonly referred to as apparent polar wander paths; APWPs) have arcing trajectories similar to fracture zones and hotspot tracks. This similar geometry is expected given that they are also tracking plate motion (Figure 1). Building on the approach of Francheteau and Sclater (1969), Gordon et al. (1984) developed an algorithm to find the best-fit small circles to paleomagnetic poles, which would furnish Euler poles for the plate in question for that time period. This approach for constraining APWPs, called paleomagnetic Euler pole analysis, has the attractive feature of providing a complete description of the plate motion, including paleolongitudinal changes and therefore rates of motion. However, it has the drawback of not providing readily computed uncertainties associated with the paleomagnetic poles and not incorporating age uncertainties. An additional complexity is that while relative motions between tectonic plates can be consistent over extended periods and explained by a single Euler, absolute plate motion relative to the spin axis can be more time variable due to changes associated with one of multiple Euler rotations through a plate circuit or due to a combination of plate tectonic motion and true polar wander. With some notable exceptions (e.g., Beck, 1989;Beck & Housen, 2003;Bryan & Gordon, 1986;Smirnov & Tarduno, 2010;Tarling & Abdeldayem, 1996), paleomagnetic Euler pole analysis has not seen wide adoption.
In this contribution, we extend paleomagnetic Euler pole analysis by placing it within a Bayesian statistical framework and demonstrate how to invert for paleomagnetic Euler poles using Markov chain Monte Carlo (MCMC) methods. This framework has the advantage of naturally incorporating uncertainties in paleomagnetic pole positions, as well as widely disparate age uncertainties associated with individual paleomagnetic poles. The resulting stage poles from these inversions are not a single answer, but are instead a distribution of possible answers, furnishing uncertainties as part of the solution process. Proper treatment of uncertainties also allow researchers to avoid overfitting models to noise. For instance, Iaffaldano et al. (2012) employed a similar Bayesian approach to inversions for finite plate rotations. They used seafloor data to reconstruct relative plate motion across Atlantic, Indian and South Pacific ridges incorporating uncertainties into the inversion. In the process, they demonstrated that interpretations of changes in relative plate motions on timescales <1 Myr are the result of overfitting the data and concluded there are fewer kinematic changes and longer stability in plate motions than suggested by some previous treatments.
In the sections below, we first review different approaches that have been applied to synthesize and interpret APWPs. We then describe the formalism of Bayesian inversions and MCMC methods that we will apply. We then describe the statistical model which we will be inverting and demonstrate the inversion of several synthetic data sets. Finally, we show case studies where we apply the method to paleomagnetic poles in order to develop estimates of paleomagnetic Euler poles and associated plate motion.

Interpretation of Apparent Polar Wander Paths
A sequence of paleomagnetic poles from the same continental block can be synthesized into an APWP, which can then be used to develop plate tectonic reconstructions and models of plate motion with respect to Earth's spin axis through time. Interpretation of these paths becomes difficult in the case of limited, uncertain, or conflicting data, and when the ages of paleomagnetic poles are poorly known. A number of approaches to developing APWPs have been developed, which we briefly review here.

Latitudinal Drift
A paleomagnetic pole provides constraints on the paleolatitude and orientation of a lithospheric block. However, due to the rotational symmetry of Earth's time-averaged geocentric axial dipole magnetic field, paleomagnetic poles do not directly constrain absolute paleolongitude (Butler, 1992). The simplest analysis of an APWP is to compare the paleolatitudes implied by successive poles for a point on a respective block. The difference in paleolatitudes gives a minimum angular distance over which the block has traveled. When this distance is compared to the age difference between the poles, such a comparison establishes a rate of latitudinal motion.
It is possible to estimate confidence bounds on the rate of latitudinal drift through bootstrap resampling (e.g., Tarduno et al., 1990) or by taking a Monte Carlo approach. Swanson-Hysell et al. (2014) developed a Monte Carlo sampling method and applied the method to a pair of poles from the Proterozoic Midcontinent Rift of North America to estimate the range of implied latitudinal drift. They also sampled from the uncertainties of radiometric dates associated with the poles, assuming Gaussian distributions, in order to incorporate age uncertainties into the analysis. With samples of pole position and ages, they were able to estimate the 95% confidence estimates on the rate of latitudinal drift.
Whether using point estimates of the latitudinal drift rate or using Monte Carlo estimates, the latitudinal drift interpretation of APWPs remains limited as it represents a minimum estimate of total plate motion. It does not resolve longitudinal drift rate, nor does it naturally extend to APWPs with more than two poles, especially if two coeval poles are not in agreement, as it requires the selection of pole pairs.

Spherical Splines
When considering APWPs with many poles, it becomes more difficult to perform latitudinal comparisons between pairs of poles. It is not always clear which pairs of poles to compare in cases where there are many overlapping paleomagnetic poles that have variable uncertainties associated with their positions and ages.
One approach to synthesize such data is to fit a spline through the set of paleomagnetic poles, constraining the path to lie on the surface of a sphere. This approach was pioneered by Torsvik et al. (1992) using the spherical spline algorithm developed by Jupp and Kent (1987). This approach has the advantage of allowing the paleomagnetic poles to be weighted by their spatial uncertainties. The uncertainty assigned to a paleomagnetic pole can be the 95% confidence interval on the pole position, but it can also be augmented by various quality screening factors, such as the quality ("Q") factor of Van der Voo (1990) (Torsvik et al., 1992). Even with the weighting of the paleomagnetic poles by uncertainty, there can be unrealistic loops in the APWP generated by the spline fit. To combat this behavior, the spline can also be computed under tension which penalizes curvature and produces a smoother path (Torsvik et al., 1996).
The spherical spline approach to interpreting APWPs has attractive features. It produces a smooth path through the data that can incorporate spatial uncertainties, and may be efficiently computed. However, it does have some drawbacks. It is not easy to determine the appropriate uncertainty weighting and spline tension parameters for the fit, and what A finite rotation of a plate around an Euler pole (dark red circle) results in arcuate oceanic fracture zones and hotspot tracks (cartoon mountains) which describe small circles on the globe. The same finite rotation produces a circle in the apparent polar wander paths (APWP), which is illustrated by blue paleomagnetic poles. By fitting a small/great circle to the APWP, we may recover the Euler pole that produced the rotation which is termed the paleomagnetic Euler pole. Cartoon is adapted from Gordon et al. (1984 (Torsvik et al., 1996). It also does not have a simple way of incorporating age uncertainties of the paleomagnetic poles. Finally, by their nature, splines do not readily represent the sharp hairpin cusps that characterize abrupt shifts in motion that plates sometimes undergo (Gordon et al., 1984;Irving & Park, 1972;Torsvik et al., 2012).

Running Means
An alternative method for developing APWPs is to perform a running Fisher mean on the poles with a moving window through time (Irving, 1977;Torsvik et al., 2008;Van der Voo & Torsvik, 2001). In such an analysis, paleomagnetic poles in a compilation are averaged with a defined window duration (typically 10-30 Myr) that controls the amount of smoothing. Like spherical splines, the running mean approach has the ability to effectively damp the effect of outlier poles that could lead to spurious motion in the APWP if there are sufficient data. Torsvik et al. (2008) also investigated the effects of combining running means with spherical splines, by first computing a set of mean poles and then fitting a spline through those means. The simplicity of using the running mean approach to develop an APWP, as well as its ability to suppress potentially spurious poles, has led to it being widely adopted.
While it is a straightforward approach with advantages for synthesizing paleomagnetic poles, the running mean approach shares many of the drawbacks of the spline approach. As with the spherical spline method, age uncertainties on the poles are not incorporated, nor are the uncertainties in pole positions as typically implemented. It is not obvious how to best choose the window duration, and different window durations are likely appropriate for different data sets. It is also unclear how to interpret the resulting uncertainties in the path that are reported as a Fisher A 95 ellipse of the mean of the poles as the poles need not be coming from a Fisher distribution particularly if there is ongoing plate motion. A recent study discussed the advantage that could be gained by calculating APWPs from site level data (virtual geomagnetic poles) rather than from study level mean paleomagnetic poles (Vaes et al., 2022). This work highlighted how such an approach could ameliorate the disadvantage of typical running mean approaches of giving equal weight to study level poles which themselves comprise varying numbers of site level virtual geomagnetic poles.

Paleomagnetic Euler Poles
The idea of fitting small circles to paleomagnetic data to gain insight into plate motions was first described by Francheteau and Sclater (1969). The center of such a circle on a sphere is the pole of rotation which was termed a paleomagnetic Euler pole by Gordon et al. (1984) who developed a method to invert for the position of such poles. The approach builds on the recognition that plate motions can often be approximated by finite rotations around Euler poles which can potentially be steady for millions or tens of millions of years. As a result, an arcuate APWP of a plate can be described by Euler rotations, which produce small circle paths on Earth's surface if the Euler pole remains in a relative constant position ( Figure 1). Since the analysis solves for the Euler pole that produces a given small circle, this approach allows for an estimate of the full motion of a given plate, including the total plate speed (instead of just the latitudinal component of the speed). Effectively, the analysis utilizes both the change in plate orientation and latitude information embedded within paleomagnetic poles to estimate both latitudinal and longitudinal motion as succinctly specified by Euler pole parameters.
However, paleomagnetic Euler pole analysis has many of the same deficiencies that spline fits and running means have-it is not easy to compute uncertainties, especially in the presence of unknown ages of poles. Furthermore, one has the additional challenge and related uncertainty of deciding how many paleomagnetic Euler poles to include for a given sequence of paleomagnetic poles. In most treatments, this number of stage poles is subject to researcher choice (however, see theoretical progress in this regard by Gallo et al. (2022)). In this contribution, we develop a Bayesian statistical approach to paleomagnetic Euler pole analysis which attempts to address some of these deficiencies.

A General Description of Inverse Problems
The central question motivating inverse problems is "How probable is a particular model, given my observations?." We represent observations by the data vector d, and a model by the vector of model parameters m, so the above question can be expressed as the function P(m|d) (the probability of the model given the data). Traditional 5 of 26 frequentist approaches to an inverse problem often proceed by maximizing the likelihood function, defined by the probability of the data given a particular model (e.g., Aster et al., 2005): (1) The likelihood function replaces something that is difficult to compute (namely, P(m|d)) with something that is less difficult to compute. To compute ( | ) we need to have two things: a statistical model for uncertainties in the observations d and a forward model that allows us to compute predictions. We denote the forward model by g: where the superscript "p" denotes a predicted value. If each of the observed data d i are described by Gaussian random variables with standard deviations σ i , the likelihood function is given by the product of the individual likelihoods of the observations: The likelihood function  is maximized by searching over the model parameter space. If the uncertainties in the observations are Gaussian, then maximizing the likelihood function is equivalent to the least squares solution (Aster et al., 2005).
A standard maximum likelihood fit will frequently overfit the observations, resulting in unrealistic solutions. In the context of APWPs, these overfit solutions may pass through every paleomagnetic pole, including less reliable ones, resulting in loopy or jerky paths. In order to address such overfitting, some form of regularization is usually included in the solution of the inverse problem, such as penalizing the magnitude or curvature of the solution. Both the running-mean and the spline under tension approaches to APWPs are forms of regularization on the problem.

Bayesian Approach
The Bayesian approach to inverse problems takes a different strategy from the frequentist one. Rather than finding a model with parameters that maximize the likelihood, it treats the underlying model as a set of random variables with individual probability distributions. The probability distribution of the model given the data (P(m|d)) is then found by an application of Bayes theorem (cf. Sivia & Skilling, 2006): It is often unnecessary to calculate the denominator of Equation 4, which is a normalization constant, leaving us with The quantity P(m|d) is known as the posterior probability, and it represents our desired knowledge about the distributions of the parameters m. The first factor on the right-hand-side of Equation 5 is identical to the likelihood function described above, and the second factor is known as the prior probability of the model. In the case for Euler pole and/or true polar wander inversions, our model parameters are Euler pole positions (longitudes, latitudes) and rotation rate as well as changepoints (when multiple rotations are inverted for).
In this case, the prior probability distributions reflect the state of our knowledge and beliefs of the values of the Euler rotation axes and rotation rates prior to the consideration of data. It allows us to incorporate theoretically or empirically derived constraints that are not otherwise included in the forward model. In contrast with the classical statistical approach of regularization, the Bayesian inverse problem can (in effect) regularize the problem by making choices of probability distributions that have less probability density in parameter space with less realistic values (e.g., Minson et al., 2013;Sambridge et al., 2013).

Markov Chain Monte Carlo Methods
For complex models, it is usually impossible to calculate the posterior probability distribution in Equation 4 directly (Davidson-Pilon, 2015). It is much more tractable to generate a Markov chain which, upon convergence, 6 of 26 generates samples from the desired posterior distribution (Gelman et al., 2013). This approach defines a class of methods known as MCMC methods.
In comparison with drawing Monte Carlo resamples from random variables independently, MCMC resampling generates a sequence of samples where the current value for each variable is dependent on its prior value. This approach allows for more efficient sampling of probability distributions with sampling that converges toward higher posterior probabilities. For an introduction to MCMC methods, the interested reader could refer to Gelman and Rubin (1996), Sambridge et al. (2013), and Davidson-Pilon (2015). A number of high-quality open source software packages for implementing MCMC models exist, including WinBUGS (Lunn et al., 2000), PyMC (Salvatier et al., 2016), and Stan (Carpenter et al., 2017). We make extensive use of PyMC in this work.

Distributions on a Sphere
In order to proceed with a Bayesian description of the problem, every parameter in the model needs to be described by some statistical distribution that determines the probability that the parameter takes a specific value. Parameters like pole ages can be described by 1D probability distributions (such as uniform or normal distributions), whereas Euler pole locations are described by 2D distributions of directional data on the surface of a sphere. We review several of these distributions here. For a comprehensive discussion of spherical probability distributions, see Fisher et al. (1987). Plots of the following distributions (uniform, Fisher, and Watson), as well as samples drawn from them, are shown in Figure 2.

Uniform Distribution
The simplest probability distribution on a sphere is the spherical uniform distribution. It has a probability density given by where ρ U is the probability density, ϕ is the longitude, and ψ is the latitude (we refer to the Cartesian unit vector x as a concise representation of ϕ and ψ). Non-uniform distributions on a sphere reduce to the uniform distribution in some limit (i.e., the Fisher distribution as the precision parameter goes to zero). We use the uniform distribution when we want to specify an uninformative prior distribution for directional parameters.

Fisher Distribution
The Fisher distribution (also called the von Mises-Fisher distribution) is the analogue of a 2D normal distribution on a sphere ( Figure 2). The probability density ρ F at a point x is given by where κ F is the concentration of the distribution, the unit vector for the mean direction of the distribution, and C F is a normalization coefficient. It can be alternatively parameterized using θ, which is the angle between x and . The normalization factor is given by When κ F goes to zero, the Fisher distribution is equivalent to the spherical uniform distribution.
The uncertainty ellipses for paleomagnetic poles are typically calculated assuming a Fisher distribution of the underlying data, and we will use this distribution to calculate the likelihood function for pole positions in the model.

Watson Girdle Distribution
Whereas the Fisher distribution concentrates probability density around a pole on the surface of the sphere, the Watson girdle probability distribution is concentrated in a belt orthogonal to the pole ( Figure 2). It is useful for characterizing planar data, and is given by where κ W is the concentration of the girdle, C W is a normalization coefficient, and the other parameters are identical to those in the Fisher distribution. The Watson distribution is girdle-shaped only when κ W is a negative number, which is the only case we consider here.
The normalization factor is given by where 1 F 1 () is Kummer's confluent hypergeometric function, which is available in most software libraries of special mathematical functions. As with the Fisher distribution, when κ W goes to zero, the Watson distribution is equivalent to the spherical uniform distribution.

Forward Model
A forward model describes how we generate predicted observations given a set of model parameters (Equation 2). The forward model for paleomagnetic Euler pole analysis in this study is essentially unchanged from that of Gordon et al. (1984). We will consider three overall scenarios of plate motions (and hence paleomagnetic pole motions). The first is that plate motion is the result of plate tectonic motion about an Euler pole or a series of Euler poles. Each Euler pole has three parameters: a latitude, a longitude, and a rotation rate. In a model with multiple Euler poles, we also must specify the ages where one Euler pole switches to the next as an additional unknown variable. In the context of parameter inversion, these ages are often known as "changepoints." For plate tectonic motions, the motion can be along any plane intersecting the sphere such that they are small circles. The second scenario we will consider is one wherein the sole cause of the motion of paleomagnetic poles is true polar wander in which case motion is restricted to be along a great circle path coaxial to the axis of minimum inertia (Creveling et al., 2012). In such a model, three parameters exist: a latitude and a longitude of the true polar wander rotation axis and a rotation rate. A third scenario is that plate tectonic and true polar wander rotations are occurring concurrently leading to an observed paleomagnetic pole path.
We also need to initiate a model with a starting position on the globe, which, in practice, can be sampled from the Fisher distribution of the oldest paleomagnetic pole in the data set. The starting point contributes two parameters (a latitude and a longitude). Therefore, a model with n e Euler rotations will have 3n e parameters for the poles, (n e − 1) parameters for the changepoints, and two parameters for the starting location. The number parameters for which we are inverting is then given by Two more parameters will be included in a model where true polar wander is considered to happen on top of n e Euler rotations. In this case, the axis of true polar wander is constrained to be on a great circle that is orthogonal to the vector representing the starting pole position. Therefore, a total of 4n e +3 parameters are inverted for.
With the axes and rates for such spherical rotations defined, one can calculate the linear velocity v of a point p on the surface of the globe associated with angular velocity ω i by Finite rotations can be performed by constructing Euler angle rotation matrices (cf. Goldstein, 1965). We generate synthetic paleomagnetic pole positions from the forward model by stringing together finite rotations through the stage poles until the age of the observed paleomagnetic pole is reached. These positions can then be compared to the observed paleomagnetic poles in a given dataset.

Choice of Prior Distributions
Bayesian analysis requires us to specify prior probability distributions for each of the model parameters in the inverse problem. These distributions reflect our state of knowledge about the values of the parameters before we begin, and allow us the option of incorporating information otherwise not captured by the model. To avoid biasing the results of the model toward a specific posterior distribution, we usually try to choose prior distributions that are as uninformative as possible. Depending upon the context, and the type of parameter, that choice may vary. The central parameters in the paleomagnetic Euler pole problem are the Euler pole positions, the Euler pole magnitudes, the changepoints, the starting point, and the paleomagnetic pole ages, which we treat in turn. We use the notation x ∼ y to indicate that the parameter x is drawn from distribution y.

Euler Rotation Vectors
The first parameter we consider is the position of the Euler poles, which should be drawn from a spherical probability distribution. The least informative prior distribution for an Euler pole position is the uniform spherical distribution:̂∼ .
essentially allowing the Euler pole to be anywhere on the globe with equal probability.
An alternative choice is to inform our prior distribution for Euler pole positions based on observation of modern plate motions. It has long been observed that, to first order, plate motions are well explained by slab-pull torques acting along subduction zones, and to a lesser extent, ridge push and mantle traction effects (Forsyth & Uyeda, 1975;Gordon et al., 1978;Richardson, 1992).
We can ask the question of whether the Euler pole for a given plate is more likely to be on top of the plate (corresponding to a spinning motion for that plate) or away from that plate (corresponding to motion across the surface of Earth). Given that tectonic plates can broadly be considered to be the surface expression of mantle convection, we can hypothesize that the second possibility is more likely because a spinning plate has no divergence (i.e., spreading centers and subduction zones, (Forte & Peltier, 1987;Gable et al., 1991)). Without divergence, the plate motion does not contribute to large-scale convection.
To evaluate this hypothesis, we generated position samples on the surface of Earth and computed the angular distance between that point and the Euler pole for the plate in which that point resides. We used the NNR-MORVEL56 model for current plate motions (Argus et al., 2011) and restricted our analysis to the 14 largest plates. We then fit those angular distance samples to a Watson girdle distribution (Equation 9), inverting for the concentration parameter κ W . If an Euler pole position has no preference for being a particular angular distance from a point on a plate, then κ W should be close to zero, corresponding to a uniform distribution. We find that the distribution is best fit with κ W ≈ −0.7, which corresponds to the Euler pole probability density being roughly twice as large 90° away from a given point than on top of the point (Figure 3). 9 of 26

Euler Rotation Magnitudes
The magnitude of each Euler pole rotation is a positive number, specifying the rotation rate about the Euler pole (negative magnitudes can be accommodated by flipping an Euler pole to the antipode). There are several possibilities for the prior distribution for the rates. In order to not bias the inversion toward a particular rate, we can choose a uniform prior distribution with large support: where U(⋅, ⋅) is a uniform distribution between two values, and is specified in degrees per million years. Typical rotation rates for present day plate motions are under 1°/Myr (Argus et al., 2011), which corresponds to plate speeds of about 11 cm/yr for a rotation about an Euler pole that is 90° from a plate.
Another option is to choose a weakly informative prior distribution for the Euler pole magnitudes informed by recent plate motions (similar our approach of the Watson girdle prior distribution for Euler pole position). Zahirovic et al. (2015) found, based on analysis of Cenozoic and Mesozoic plate reconstructions, that plate speeds much higher than 15 cm/yr were unlikely to be sustained. A reasonable choice of distribution for strictly positive numbers is the exponential distribution, given by which has higher probability density at lower values, and falls off exponentially toward higher values. We sampled the current plate rates on Earth's surface according to NNR-MORVEL56 and fit those to an exponential distribution. The best fitting scale parameter λ for current plate rates is ∼2.5 ( Figure 3). Making this choice of prior distribution for Euler pole rotation rates can be seen as a form of regularization on plate speeds. A smaller for λ (such as λ = 1 as is shown in Figure 3) could reflect similar knowledge about the distribution of plate speeds with a less restrictive regularization. The histogram is the angular rotation rate from one thousand samples from the surface of Earth, using the NNR-MORVEL56 model. A fit to this sample set with an exponential distribution yields a scale parameter of λ ≈ 2.5. We also show the distribution for λ = 1.0, which imposes less regularization on the rate as a less restrictive prior distribution, and a uniform U(0, 4) distribution between 0 and 4°/Myr, which specifies no preference for slower speeds if set as the prior probability. (b) Prior probability density for the position of the Euler poles, with the north pole as the site latitude and longitude. We again sampled one thousand points on Earth's surface, calculating the angular distance between that point and the Euler pole for its plate. If we model the probability distribution as being drawn from a Watson distribution, these angular distances correspond to colatitudes, where the pole is the sampled point. Fitting the resulting angular distribution to a Watson girdle distribution finds κ W ≈ −0.7. Since the Watson distribution is rotationally symmetric, longitudes do not contribute to the fit. For κ W ≈ −0.7 the probability density is roughly twice as large at the equator (90° from a plate) as at the pole (on top of the plate).

Changepoints
Changepoints occur sequentially between the oldest (at age a max ) and youngest (at age a min ) paleomagnetic poles.
We choose a uniform distribution as a prior for these changepoints: where c i is the i'th changepoint.

Starting Position
Finally, the starting position xstart for the set of Euler pole rotations needs a prior distribution. We could choose another uniform distribution, but a more reasonable choice is to start near the oldest paleomagnetic pole in the dataset. We therefore choose the Fisher distribution of the oldest paleomagnetic pole as a reasonable prior distribution for a start point:x where κ F0 and 0 are the concentration parameter and mean direction of the oldest paleomagnetic pole in the dataset.

Pole Ages
One of the major advantages of Bayesian analysis is the ability to naturally incorporate uncertainties in as many parameters as needed. Previous approaches to modeling APWPs have the drawback that they do not easily account for uncertainties in the age of paleomagnetic poles. In our approach, we can include age uncertainty by including the age of the poles and associated uncertainty as parameters in our model.
There are many different ways to constrain the ages of the geologic units from which we obtain paleomagnetic poles, including radiometric dating, biostratigraphy, magnetostratigraphy, and cross-cutting relationships. Here we concentrate on poles that are either interpreted to be the age of a single radiometric date or are interpreted to be bracketed stratigraphically between two dates (derived radiometrically or by using other age control such as biostratigraphy). If a geologic unit has been radiometrically dated, we can model the age of the j'th paleomagnetic pole a j as a normal distribution with mean μ j and standard deviation σ j : where N(., .) denotes a normal distribution.
Frequently, however, the geologic unit from which we obtain a paleomagnetic pole is not well dated, but its age can be constrained to lie between those of well-dated units stratigraphically above and below it, dates obtained by cross-cutting relationships, or a number of dated units within a broader province. In these cases, a uniform distribution between those ages is a reasonable choice for the prior distribution: where a young and a old are the ages of the lower and upper age constraints, respectively.

True Polar Wander
Another advantage of this method is the ability to also incorporate a component of true polar wander rotation on top of Euler pole rotations. Because true polar wander represents the rotation of the entire solid Earth with respect to the spin axis, an APWP that is solely consisted of true polar wander should follow a great circle trajectory and the associated rotation axis should be orthogonal to the great circle. In a case where both small circle rotations and true polar wander are recorded by an APWP, the true polar wander rotation axis is constrained to be 90° away from the start position of the path and the rotation rate can be defined in a similar fashion as that of the Euler pole magnitudes.
To summarize our choices for prior distributions: • Euler pole positions: spherical uniform distribution, or a Watson girdle distribution with κ W ≈ −0.7 or a value closer to 0.
ROSE ET AL.
10.1029/2021JB023890 11 of 26 • Euler pole magnitudes: Uniform distribution, or an exponential distribution with λ ≈ 2.5 or smaller. • Changepoints: uniform distribution between a min and a max . • Starting position: Fisher distribution usually defined by the Fisher statistics of the oldest paleomagnetic pole position in an APWP. • Paleomagnetic pole ages: normal or uniform distribution, depending on the type of age control for the geologic unit from which the pole was obtained.

Likelihood
In addition to the choice of prior distributions, we need a statistical description of the observations. This description will allow us to calculate the likelihood function, which, when combined with the prior distributions, allows us to evaluate Bayes' theorem (Equation 5).
In the case of APWPs, our observations are paleomagnetic poles. The most common statistical distribution for describing paleomagnetic poles is the Fisher distribution (although other distributions are sometimes used, such as the Kent or Bingham distributions, c f. Tauxe (2010)). Given the set of model parameters m and the forward model g(m), described above, we can calculate the predicted paleomagnetic pole unit vectors x . For a set of n paleomagnetic poles, the likelihood is then given by the product of the probabilities of each observation of paleomagnetic pole position:

Example Inversions
Before proceeding with inversions for paleomagnetic Euler poles using real paleomagnetic data, it is useful to consider a few examples of inversions for idealized synthetic datasets. We have specified the forward model described above using the package PyMC (Salvatier et al., 2016) which enables us to perform the inversion. Within PyMC, we are able to specify custom probability distributions enabling the directional data distributions illustrated in Figure 2 to be implemented. The MCMC analysis uses the Metropolis-Hastings algorithm for sampling as it can be applied to custom probability distributions in contrast with other sampling algorithms within PyMC. Our code for the inversions has an open-source GPL license and is available on Github (https://github.com/Swanson-Hysell-Group/Bayesian_PEP_inversion) and archived on Zenodo (https://doi.org/10.5281/zenodo.7087845).

One Euler Rotation
We begin by trying to recover the Euler pole for a single rotation. We generate an idealized synthetic APWP of four poles by starting from a pole at 19°N, 024°E, and rotating around an Euler pole at 00°N, 000°E for 100 Myr at a rate of 1°/Myr. We produce paleomagnetic poles at 100 , 75 , 50 , and 25 Ma, and prescribe A 95 of 4° to each pole (where A 95 indicates the 95% angular confidence for the pole position).
To demonstrate the ability to recover the Euler rotation axes and rotation rate using the inversion method, we introduce minimal prior knowledge on these parameters, but assume the ages of the paleomagnetic poles are well constrained. Therefore, we use a Watson girdle distribution with κ W = −0.1 (more weakly informative than −0.7 of Figure 3 and approaching a uniform distribution as illustrated in Figure S2 in Supporting Information S1) as the prior distribution for the Euler rotation axis, and a uniform distribution of U(0, 4) as the prior distribution for the Euler rotation rate (in °/Myr).
The results of the inversion are shown in Figure 4. The Bayesian approach successfully recovers a posterior probability distribution for the position of the Euler pole, as well as a rate that is centered near the true value of 1°/Myr (Figure 4b). The posterior distribution for the rate has a highest posterior density credible interval at 95% (which we abbreviate from here as a 95% credible interval) between 0.7°/Myr and 1.2°/Myr, reflecting the resolving power of the inversion for data with the given uncertainties. An example of this same inversion with simulated noise on the age and position of paleomagnetic poles is developed in PEP_synthetic.ipynb within the archived code repository and shown in Figure S1 in Supporting Information S1. It is also successful in recovering the true values within the credible intervals of the posterior distributions.

Two Euler Rotations
We next consider an inversion for an APWP with two stage poles. Unlike the previous example where we inverted for a single Euler pole, this inversion also requires a changepoint. We generate nine idealized paleomagnetic poles with A 95 of 4° for each from a starting point at 5°S, 030°E. The first rotation is around an Euler pole at 10°S, 000°E, and rotates at 1.5°/Myr for 60 Myr. The second rotation is around an Euler pole at 10°N, 300°E and rotates at a rate of 0.75°/Myr for the same amount of time. The synthetic poles associated with this two stage model are shown in Figure 4.
During the inversion, the prior knowledge for the Euler rotation axes and rates and the age of the paleomagnetic poles are set in the same way as in the one Euler rotation section above. The prior distribution for the changepoint is set as a uniform distribution with the minimum and maximum ages being the age of the youngest and oldest paleomagnetic poles. The inversion successfully recovers the Euler pole rotation rates with posterior distributions centered near the true values. The inversion also successfully recovers the changepoint of 60 Ma between the first and second Euler pole with a 95% credible interval for the changepoint of 65 to 55 Ma. The inverted Euler pole positions for the first stage pole (blue distributions in Figure 4c) are centered near the true location (blue star in Figure 4c) and the inverted second Euler pole positions (red distributions in Figure 4c) encompass the true location as well (red star in Figure 4c). There is a larger spread in credible Euler pole positions for the second Euler pole particularly in the direction perpendicular to the APWP. The greater uncertainty on the location of the second Euler pole results from the motion along the APWP being smaller relative to the uncertainty on the paleomagnetic pole positions. As a result, the inverted paths have more variable curvature.
Overall, this example demonstrates that the inversion framework can resolve APWPs with more than one stage pole and provides posterior probability distributions on the Euler rotation changepoint in addition to the multiple Euler poles positions and rates (Figures 4c and 4d).

Incorporating Age Uncertainty
A benefit of the Bayesian approach to inverse problems is its generality. As long as some effect can be described statistically and incorporated into our forward model, we can include it in the inverse problem. As a result, we are able to include uncertainties on the ages of paleomagnetic poles. Assigning a fixed single age to paleomagnetic poles can bias an analysis when the age of a pole is not precisely known. The inversion framework enables the varying certainty on pole ages to be incorporated into the analysis. To demonstrate this capability, we use a similar test case as in the one Euler pole inversion with the ages of synthetic poles being adjusted to span from 140 to 40 Ma (poles of 140, 115, 90, 65, and 40 Ma) (Figure 5a), but assign prior distributions for the ages of the poles. For the first and last poles, we assume they are radiometrically dated with standard deviations of 5 Myr. However, we assume that the middle three poles have no age control, except that their ages are constrained to be between the first and last poles. We thus assign Gaussian prior distributions to the first and last poles and uniform prior distributions to the middle three (Figure 5b). Despite this uninformative prior distribution on the age of the middle three poles, the inversion successfully places the estimated distributions of ages of the middle three poles to be centered at ca. 115 , ca. 90 , and ca. 65 Ma as can be seen in their posterior age distributions (Figure 5d).
For real data, adding uncertainties to the ages of the poles enables us to properly represent our knowledge of the constraints on the APWP. These uncertainties enable data to constrain the location of the path without providing an overly tight constraint on the timing when the true age is uncertain. The resulting posterior distributions provide predicted ages for the poles associated with Euler pole inversions (Figure 5d).

Reporting the Apparent Polar Wander Path
A difficulty with the Bayesian approach is that the credible interval of Euler pole posterior positions are not easy to report as they do not neatly correspond to a parametric statistical distribution. We visualize the solutions with spatial histograms for the Euler pole positions and by showing example inverted paths. A typical product that one is seeking with such an inversion is an apparent polar wander path reported as interpreted pole position at given intervals. An approach that can be taken is to calculate a number of predicted pole positions at a given time implied by inverted Euler pole models and then to calculate the Fisher mean of these inverted pole positions as was done in Swanson-Hysell et al. (2019). This approach provides the mean pole position on the apparent polar wander path. We also wish to report the uncertainty on that estimated position. The spread in the position of these inverted pole positions has real meaning related to the certainty of the path position at a given time resulting from the Euler pole positions and rotation rates in the posterior distribution. The spread in pole positions can be approximated as a Fisher distribution from which the angle from the mean that contains 95% of the solutions (θ 95 ) can be calculated: This angle is analogous to a 2σ uncertainty in Gaussian statistics. The pole positions typically are consistent with being drawn from a Fisher distribution such that reporting the 95% confidence of angular deviation can be appropriate, but note that the distributions are not necessarily Fisherian such that this can be a simplified approximation. Regardless, implied pole positions from these inversions are more tightly clustered than the inverted ROSE ET AL.

10.1029/2021JB023890
14 of 26 Euler pole positions themselves. This approach for summarizing an APWP is applied to Australia's Cenozoic path in the following section.

Application to Australia's Cenozoic APWP
Australia and East Antarctica were connected from the Proterozoic Eon until the Mesozoic Era when rifting led to the onset of seafloor spreading between the continents ca. 83 Ma (Veevers, 2012;Williams et al., 2011). Australia's plate motion relative to East Antarctica is considered to be well constrained from ca. 61 Ma to the present on the basis of fracture zones and magnetic anomaly data (Cande & Stock, 2004). The plate motion model of Müller et al. (2016) utilizes the reconstruction of Australia relative to East Antarctica of Cande and Stock (2004) for 38 Ma to the present-day and that of Whittaker et al. (2007) for 100 to 38 Ma. The absolute motion of Australia is constrained in the model through a plate circuit that reconstructs Australia relative to East Antarctica, East Antarctica relative to Africa, and Africa relative to the spin axis (Müller et al., 2016). The relative motions between the plates in this plate circuit that are implemented in the plate model of Torsvik and Cocks (2017) have small differences from that of Müller et al. (2016) such as the use of Tikku and Cande (2000) for the Australia to East Antarctica rotation between 84 and 44 Ma. Despite these small differences, the overall relative motions are quite similar. A more significant difference is the use of a different reference frame for the solid Earth relative to 15 of 26 spin axis between Müller et al. (2016) and Torsvik and Cocks (2017). The choice in this regard is important as pertains to comparisons with paleomagnetic data as paleomagnetic data is in the reference frame of the spin axis. The model of Müller et al. (2016) reconstructs plates relative to the global moving hotspot reference frame of Torsvik et al. (2008) which is taken as fixed relative to the spin axis over the Cenozoic. In contrast, Torsvik and Cocks (2017) utilize the reference frame of Doubrovine et al. (2012) which includes true polar wander rotation between the global moving hotspot reference frame and the spin axis. Notably, a recent analysis of interpreted hotspot tracks in Eastern Australia by Hansma and Tohver (2020) argued that their position in a paleomagnetic reference frame is consistent with the true polar wander implemented in the Doubrovine et al. (2012) reference frame. In Figure 6, we show the implied pole position for Australia that results from these global plate motion models. Both models imply faster Euler rotation rates after ∼35 Ma than before that time associated with an acceleration in Australia's northward drift (Figure 6d). We use these plate motion models that are largely based on seafloor data as a comparison to an inversion of Cenozoic paleomagnetic data from Australia for the past 60 million years.
There is a contentious history of interpreting Australia's Cenozoic paleomagnetic record leading to varying APWPs with discussions of the relative fidelity of igneous and sedimentary poles in the literature (e.g., Hansma & Tohver, 2019;Idnurm, 1985;Idnurm, 1994;Musgrave, 1989). The paleomagnetic pole database for Australia in the Cenozoic has improved substantially with the development of both new Oligocene and Miocene paleomagnetic data (Hansma & Tohver, 2018 and Ar-Ar dates that provide ages for both previously undated units and supersede previous K-Ar dates (Cohen et al., , 2013Knesel et al., 2008). This updated paleomagnetic pole database is listed in Table 1 and visualized in Figure 6a. These paleomagnetic poles give us the opportunity to develop paleomagnetic Euler pole inversions that can be compared to the independent plate tectonic reconstructions of Müller et al. (2016) and Torsvik and Cocks (2017).
The paleomagnetic data necessitate at least one change in the Euler pole between 60 Ma and the present. This need for more than one Euler pole is most readily visualized in the latitude of the poles which implies a change in rate between 60 Ma and the present (Figure 6b). We apply the Bayesian inversion framework to invert for two Euler poles for Cenozoic Australia using prior probability constraints from these paleomagnetic poles. We use an uninformed prior probability for the timing of the changepoint between the two Euler poles in which it is assigned a uniform probability distribution between the age of the oldest (∼60 Ma) and the youngest paleomagnetic poles (∼0.005 Ma) in the compilation. We also use a weakly informative uniform prior distribution of U(0, 4) for the Euler rotation rate, as we know based on the position and age of the paleomagnetic poles that the overall plate Note. PLat = pole latitude; PLon = pole longitude; A 95 = 95% angular confidence bounds on pole positions; N = number of site level VGPs used to calculate the mean pole positions; Pmag reference = reference for paleomagnetic pole data; Age (Ma) = nominal age calculated from age constraints on paleomagnetic poles; Age min = lower bounds for pole ages; Age max = upper bounds for pole ages; Age information = summary description for geochronology methods for pole ages and associated references; Dist type = choice of using uniform or normal distributions for assigning prior uncertainties for pole ages. motion rate for Australia in the Cenozoic is most likely to be less than 4°/Myr (Figure 6a). For the prior distribution for the Euler rotation axes, we use a weakly informative Watson girdle distribution with κ W = −0.1.
The posterior distribution of the Euler pole positions, samples of the small circle paths generated from the posterior distributions, the implied pole latitudes, and the full plate motion angular rotation rates are shown in Figure 6. The probability density of the inverted positions for the younger Euler pole (ca. 40 Ma to present) is more tightly constrained than the older Euler pole position given both the higher number of paleomagnetic poles and longer track length. The posterior distribution of the younger Euler pole positions (red density in Figure 6b) recovered through the inversion correspond well with the Euler pole positions extracted from the Müller et al. (2016) and Torsvik and Cocks (2017) models (red and pink stars in Figure 6c). The inversion is also successful in recovering the change in rate seen in both models where Australia's plate motion accelerated in the later portion of the Eocene (Figure 6b). The posterior distributions for Euler rotation rates have median values that increase from 0.3°/Myr to 0.8°/Myr (Figure 6d). The 95% credible interval for the age of the changepoint associated with this change in rate spans from ca. 51 to 30 Ma with a median 40.8 Ma. While this broad posterior distribution reflects the sparse paleomagnetic poles constraints between 60 and 40 Ma, the timing of the change in rate is consistent with the Müller et al. (2016) and Torsvik and Cocks (2017) models ( Figure 6b).
There is a closer correspondence between the latitude of the path inverted using the paleomagnetic pole constraints with the Müller et al. (2016) model than the Torsvik and Cocks (2017) model as can be seen by the inverted paths (red for younger Euler; blue for older Euler) plotting atop the thicker dashed line in Figure 6b. In contrast, there is a tighter correspondence between the longitude of the pole path resulting from the inversion with the Torsvik and Cocks (2017) model as seen in Figure 6c where the inverted paths are atop the squares which mark the Torsvik and Cocks (2017) path. The path of Müller et al. (2016) deviates to the east from the inverted paths whose trajectory are being constrained by the ca. 60 Ma North Rankin 1 Drillcore pole (Table 1). To facilitate these comparisons, we apply the approach described in the Reporting the apparent polar wander path section to summarize the new APWP (Table 2)  It is intriguing that the inclusion of the Doubrovine et al. (2012) true polar wander correction in the Torsvik and Cocks (2017) model results in improved correspondence in the terms of pole longitude, but more of a discrepancy with the paleomagnetic poles and their inversion in terms of latitude (Figures 6 and 7). While the conclusions that can be drawn in this regard are limited given the restriction of the analysis to Australia, such comparisons could  be a fruitful future direction to use for comparisons between reference frames. Overall, the analysis highlights a broad agreement in the interpretation of kinematics that results from inversion of the Cenozoic paleomagnetic poles for Australia and that which results from seafloor-based plate motion models.

Application to the Keweenawan Track
An impetus for the development of this Bayesian paleomagnetic Euler pole inversion method was to constrain the absolute rates of plate motion associated with the ca. 1,109 to 1,070 Ma Keweenawan Track of paleomagnetic poles from the Proterozoic continent of Laurentia (Halls & Pesonen, 1982;Swanson-Hysell et al., 2009. These poles are associated with the rapid motion of Laurentia toward the equator leading up to the Grenvillian orogeny and the assembly of the supercontinent Rodinia (Swanson-Hysell, 2021). While previous approaches had established estimates for the rate of latitudinal motion associated with the Keweenawan Track (Davis & Green, 1997;Swanson-Hysell et al., 2014), this Bayesian inversion method was implemented in Swanson-Hysell et al. (2019) in order to constrain absolute rates (without providing the level of detail on the methodology i.e., elucidated in the present contribution). In addition to our goals of constraining absolute rates, this approach was particularly appropriate for the Keweenawan Track as it includes poles with quite disparate precision on their age constraints as some are constrained tightly by high-precision U-Pb dates (e.g., Fairchild et al., 2017) while others have looser stratigraphic constraints. A technical note related to the analysis in Swanson-Hysell et al. (2019) is that it used a version of the code written in PyMC2 (Patil et al. (2010); https://github.com/ ian-r-rose/mcplates). The PyMC project has migrated to PyMC3 (now known simply as PyMC) which is a total rewrite of the module (Salvatier et al., 2016) and necessitated our code to refactored into its present PyMC version (https://github.com/Swanson-Hysell-Group/Bayesian_PEP_inversion). These two versions of the code apply the same methodology and recover very similar solutions. One change that we adopted in this study is that

of 26
we applied inclination shallowing corrections to the sedimentary paleomagnetic poles within the Keweenawan Track (i.e., the Nonesuch Formation and lower Freda Formation poles of Henry et al. (1977)) given how widespread inclination shallowing is for remanence held by detrital hematite. While more research is needed on these sedimentary rocks to quantify the degree of inclination shallowing, we take the approach of Domeier et al. (2012) in utilizing an assumed flattening factor of 0.6 given that sample-level data are not available. A ripe direction for future work would be to adopt a probabilistic framework for inclination shallowing correction rather than an assigned single value.
For all Keweenawan Track inversions, we use weakly informative Watson girdle distributions with κ W = −0.1 for the prior distribution of Euler poles, and use an exponential distribution with λ = 2.5 for the prior distribution of Euler rotation rates (based on the analysis shown in Figure 3). Paleomagnetic Euler pole inversions to the Keweenawan Track are shown in Figures 8 and 9. The posterior distribution of the Euler pole positions are shown along with a sampling of pole positions generated from the posterior Euler parameters and pole ages which are plotted over the paleomagnetic poles (as compiled in Swanson-Hysell et al. (2019)). The prior distributions for the ages of the poles are shown in Figure 9b. The prior probabilities assigned to the ages illustrate the variable age uncertainties associated with the paleomagnetic poles in the Keweenawan Track some of which are tightly constrained by radiometric dates and others which have looser age constraints (Swanson-Hysell et al., 2019). The posterior distribution of the angular rotation rates for Laurentia resulting from the inversions are shown as angular velocities along with their 95% credible intervals (Figures 8 and 9). To compare models with varying combinations of Euler poles and assess whether the results of the inversions represent good fits for the observed paleomagnetic poles, we calculate the distribution of the natural log of the likelihood for the paleomagnetic poles generated from the posterior Euler parameters and pole ages based on the observed paleomagnetic poles. A summary of the posterior likelihood distributions is shown in Figure 8b.
Applying a single paleomagnetic Euler pole inversion to the entirety of the Keweenawan Track results in a median rate of 3.1°/Myr with a 95% credible interval of 2.7-3.6°/Myr (Figure 8). In the literature, the Keweenawan Track has variably be interpreted as being the result of fast differential plate motion (e.g., Davis & Green, 1997) or as the result of true polar wander (e.g., Evans, 2003). True polar wander results in rotation of the entire silicate Earth about an Euler pole 90° from the path and has the potential to progress at rapid rates (Rose & Buffett, 2017). Overall, true polar wander is a difficult signal to disentangle from plate motions, since any given APWP can be the result of true polar wander, plate tectonics, or some combination of the two. However, given the nature of a true polar wander rotation, an Euler pole that describes true polar wander will be 90° from the spin axis such that the Euler pole for a true polar wander event will be 90° from the paleomagnetic poles. Therefore, the extent to which true polar wander could explain motion of the Keweenawan Track can be evaluated by constraining the location of the true polar wander Euler pole to be within a great circle 90° from the path. This constraint is imposed in the inversion by using a restrictive Watson girdle distribution. As shown in Swanson-Hysell et al. (2019), a single true polar wander rotation is a poor fit to the path given its curvature. The posterior true polar wander rotation parameters results in a path that deviates from the observed paleomagnetic poles, as is illustrated by the resampled poles for younger ages falling outside of the A 95 uncertainty ellipses of the observed paleomagnetic poles in Figure 8a. As a result, the posterior log likelihood values for the true polar wander model are much smaller than those for the one Euler pole model (Figure 8b).
The pole path could be the result of a combination of true polar wander and differential plate motion. This scenario can be explored through models that invert for a true polar wander Euler pole as well as differential plate tectonic Euler poles. The model with both a plate tectonic Euler pole rotation and a concurrent true polar wander rotation partitions the motion between both with significantly faster plate tectonic motion (Figure 9a). In addition, because the true polar wander involves the rotation of the solid Earth with respect to the spin axis (rotating both the plate an its associated Euler pole), its results in more uncertainty associated with the location of the plate tectonic Euler pole, broadening the posterior Euler pole distribution (Figure 9a). Although such a rotation is possible, adding a true polar wander component to either the one or two plate tectonic Euler pole inversions does not significantly improve fits to the observed paleomagnetic poles (Figure 8b). Overall, forward models involving small circle Euler rotations yield better fits to the observed paleomagnetic data than the model that considers only true polar wander great circle rotation (Figure 8). These results indicate that there was rapid plate motion at the time which is associated with subduction that led to the closure of the Unimos Ocean leading up to collisional Grenvillian orogenesis and the associated formation of the supercontinent Rodinia ( Figure 8. For the plate tectonic Euler pole, the position (in blue density) is shown relative to the start position of the path. The true polar wander rotation will progressively rotate these plate tectonic Euler poles such that the position will change through time. The log likelihood of these models are compared with others in Figure 8b. -Hysell et al., 2022). In comparison with models that involve one Euler pole rotation, the two Euler pole inversion results in an improvement in fit as can be visualized by the modeled pole positions shown in Figure 9a and by the calculated posterior likelihood distributions for the paleomagnetic pole positions ( Figure 8b). The two Euler pole inversion results in a ca. 1,097 Ma changepoint age (with a 95% credible interval from 1,098.8 to 1,095.4 Ma) with rates that slow from ∼3.1°/Myr to ∼1.8°/Myr after the changepoint. This change could be associated with initial onset of collisional Grenvillian orogenesis along Laurentia's margin (Swanson-Hysell et al., 2019).
The results of this Bayesian Euler pole inversion method can not only provide insights into geodynamic processes, but also provide geochronologic estimates on rock units that have poor or none radiometric age constraints. As our knowledge of paleomagnetic pole ages and their uncertainties are incorporated into the inversions, the more informative prior distributions that we implement for poles that have well-constrained radiometric dates (such as the normal distributions for many igneous units of the Keweenawan Track) can act as anchor points and result in more tightly constrained posterior distributions for ages of other poles that have less informative priors (such as the uniform distributions for some igneous and the sedimentary poles of the Keweenawan Track). As is illustrated in Figure 9b, the poles whose ages could only be bracketed as uniform distributions in the prior end up having posterior age distributions with more constrained credible intervals. In this example, these ages estimated through paleomagnetic Euler pole inversion model give insight into the chronology of magmatism and basin development within the North American Midcontinent Rift.

Challenges and Opportunities in the Application of Paleomagnetic Euler Pole Analysis
While the examples in this study highlight insights that can be gained through applying paleomagnetic Euler pole analysis to segments of apparent polar wander paths, challenges remain in its broad application. One set of challenges relates to the resolving power of the method. For slow rates of motion, the position of inverted Euler poles are uncertain and it is difficult to resolve the timing of changepoints with confidence. However, this issue is more related to the precision of available constraints and is the case regardless of method. Of note is that the spread between resulting paths is much less than the spread of the posterior Euler pole position (as in the older Euler pole rotation in the Australia application; Figure 6). Determining the number of change points to implement in an inversion is an additional challenge, and one that is also more difficult for slow motions or uncertain data. This issue is addressed theoretically in Gallo et al. (2022) who apply a frequentist approach to paleomagnetic Euler pole inversion and propose a graphical approach to evaluate for the optimal number of inverted Euler poles.
The extension of such an approach beyond ideal synthetic data to noisier suites of paleomagnetic data as well as integration with our Bayesian inversion approach could prove fruitful.
It could be interpreted as an advantage of paleomagnetic Euler pole analysis that it utilizes a physical model of plate motions to inform an APWP shape that fits with our knowledge of plate motions. However, uncertainty remains regarding underlying assumptions including that Euler pole positions and associated plate motions are relatively stable over millions of years. A complexity here is that while data constrain the relative motion between plates to be relatively constant (e.g., Müller et al. (2016)), the motion of a plate relative to the spin axis can be the combined effect of multiple Euler poles in a plate circuit as well as true polar wander. Furthermore, the method also assumes that the timescale over which Euler pole positions themselves migrate are relatively rapid such that they can be approximated by instantaneous changepoints. Evaluating the validity of these assumptions is an area ripe for future research. There is the potential to apply paleomagnetic Euler pole inversion itself to gain insight into these uncertainties.
An additional opportunity comes from the method giving posterior distributions for pole ages. There is a long history of paleomagnetists comparing undated (or poorly dated) paleomagnetic poles with APWPs to gain insight on their age (Hnatyshin & Kravchinsky, 2014;McCabe et al., 1984). These methods typically rely on comparison between designated pole pairs or consider the reference APWP to not have uncertainty. In the Bayesian paleomagnetic Euler pole inversion method developed here, undated poles ( Figure 5) or poorly dated poles (Figure 9) can be incorporated into the analysis with the resulting posterior distribution age distribution providing an estimate of the age of magnetization including uncertainty.
ROSE ET AL.

Conclusions
We have extended the paleomagnetic Euler pole analysis of Gordon et al. (1984) by placing it within a Bayesian framework. As a Bayesian inverse problem, MCMC numerical methods can be used to obtain estimates of Euler pole positions and rates as constrained by paleomagnetic poles. Regularization of the inversions is not accomplished by smoothing parameters, but can instead be accomplished through prior probability distributions for the Euler pole parameters, which have clear physical interpretations. The approach enables uncertainties on both the positions and ages of paleomagnetic poles to be incorporated into the analysis. Multiple Euler poles can be included with the timing of changepoints between them being solved as part of the inversion. An advantage of this approach is that the paleomagnetic Euler poles provide an estimate for the total plate velocity including both latitudinal and relative longitudinal motion. The resulting posterior distributions from the inversions provide uncertainties for the model parameters-including estimates of plate velocity.

Data Availability Statement
The code implementing the analysis and developing the visualizations of this study are documented within Jupyter notebooks (Kluyver et al., 2016) that are available on Github (https://github.com/Swanson-Hysell-Group/ Bayesian_PEP_inversion) and Zenodo (https://doi.org/10.5281/zenodo.7087844). These notebooks utilize the functions within the Bayesian_pep library within the repository and data that are within the repository as well.
The computational environment used for running these notebooks and conducting the analysis can be reproduced using the environment.yml file within the repository. Instructions associated with recreating this computational environment can be found in the README.md file within the repository. MCMC analysis was performed using PyMC (Salvatier et al., 2016). Additional analysis was performed using the Python packages NumPy (Harris et al., 2020), Scipy , and PmagPy (Tauxe et al., 2016). Figures were created using Matplotlib (Hunter, 2007) and Cartopy (Met Office, 2010-2015.