chapter2.tex

\chapter{Methodologies}\label{chap:Metho}
\textit{This chapter mainly describes the development of a new apparent polar
wander path (APWP) similarity measuring tool used in the rest of this thesis.
APWPs based on paleomagnetic data are the principal means of describing plate
motions through most of Earth history. Comparing the spatio-temporal patterns
and trends of APWPs between different tectonic plates is important for testing
proposed paleogeographic reconstructions of past supercontinents. However, thus
far there is no clearly defined quantitative approach to determine the degree of
similarity between APWPs. This chapter proposes a new method of determining the
degree of similarity between two APWPs that combines three separate difference
metrics that assess spatial separation of coeval poles, and similarities in the
length and bearing of coeval segments, using a weighted linear summation.
Bootstrap tests are used to determine whether the differences between coeval
poles and segments are significant for the given spatial uncertainties in pole
positions. An additional Fit Quality score is used to discriminate between low
difference scores caused by comparing poorly constrained paths with large
spatial uncertainties from those caused by a close fit between well-constrained
paths. The individual and combined metrics are demonstrated using tests on
synthetic pairs of APWPs with varying degrees of spatial and geometric
difference. In a test on real paleomagnetic data, we show that these metrics can
quantify the effects of correction for inclination shallowing in sedimentary
rocks on Gondwana and Laurussia's 320\textendash0 Ma APWPs. For an APWP pair,
when one APWP's three individual metrics are all greater than or equal to, or
less than or equal to the other one's, weighting is dispensable because the
similarity ranking order becomes straightforward; otherwise assigning equal
weights is recommended, although then decision makers are allowed to arbitrarily
change weights according to their preferences.}
\vfill
\minitoc\newpage

\section{Introduction}

Paleomagnetism is an important source of information on the past motions of the
Earth's tectonic plates. The orientation of remanent magnetisations acquired by
rocks during their formation record the past position of the Earth's magnetic
poles. In older rocks, these virtual geomagnetic poles often appear to be
increasingly offset from the modern day geographic poles. Because the Earth's
geomagnetic field appears to have remained largely dipolar and centered on the
spin axis for at least the last 2 billion years~\citep{E06}, this divergence is
interpreted as recording the translation and rotation of a continent by the
motion of tectonic plates in the time since the rock formed. An Apparent Polar
Wander Path (APWP) is a time sequence of paleomagnetic poles (or, more commonly,
mean poles that average all regional paleopoles of similar age) that traces
the cumulative motion of a continental fragment relative to the Earth's spin
axis.

Investigations of the Earth's past tectonic evolution and paleogeography often
involve comparing APWPs. For example, if two now separated continental fragments
were once part of the same supercontinent, their APWPs should share the same
geometry during the interval that this supercontinent existed. If the
supercontinent has been correctly reconstructed, the APWPs should also overlap
during this interval (Fig.~\ref{fig:nambal_same_geometry}). APWP comparisons can
be used to assess plate motion models generated using different datasets and/or
fitting techniques~\cite[for example]{B02,B07,S07,T08,D11}; significant
deviations from the known APWP for a continent can also be used to identify
local tectonic rotations~\cite[for example]{G10,Ch13}. Despite the clear
importance of measuring APWP similarity, these comparisons remain largely
qualitative in nature, involving visual comparisons of specific APWP segments
and checking if they have overlapping 95\% confidence limits~\cite[for
example]{B02,B07,G10,D11}. Where quantitative measures are used, the mean great
circle distance (GCD) between coeval poles on the APWP pair has been commonly
used as a generalised difference metric for a pair of APWPs, with a lower score
indicating that they are more similar~\cite[for example]{S07,T08}. However,
because GCD is simply a measure of spatial separation and does not incorporate
geometric information about the two paths being compared, it is possible for
pairs with clearly different similarities to have similar mean GCD scores
(Fig.~\ref{fig:QualitativelyDifferentSameGCD}). Due to the inherent time
variability of the geomagnetic field, uncertainties arising from the sampling
and measurement of remanent magnetisations, and uncertainties in the
magnetization age, the mean poles that make up an APWP also have associated
spatial uncertainties. The significance of a GCD score is therefore not
immediately obvious. A score that indicates a relatively large difference
between two paths may not be significant if the spatial uncertainties are large;
a small difference could be significant if the spatial uncertainties are small
(Fig.~\ref{fig:GCDlargeButIndistinguishable}).


\begin{figure*}[tbp]
\centering
\includegraphics[width=1.01\textwidth]{../../paper/tex/ComputGeosci/figures/nambal_same_geometry.pdf}
\caption[Parts of APWPs of supercontinent fragments share the same
geometry]{(a) The APWPs for North America (black) and Baltica (grey) are
spatially distinct, but their Late Paleozoic\textendash{}Early Mesozoic sections
are geometrically similar due to them both being part of the supercontinent
Pangaea. (b) Reversing the opening of the Atlantic Ocean by rotation around a
reconstruction pole (blue star) results in the overlap of these two APWPs
between 390 million years ago (Ma) and 220 Ma, validating the proposed
paleogeography. The effects of this rotation on Baltica and its APWP (BAL) are
illustrated by the motion of the circle marker (before: blank center; after:
dark center), respectively. General Perspective projection. APWPs and rotation
parameters from~\citet{To16}.}\label{fig:nambal_same_geometry}
\end{figure*}

\begin{figure*}[!ht]
  \captionsetup[subfigure]{labelformat=empty,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{1.01\textwidth}
    \centering
    \includegraphics[width=1.01\linewidth]{../../paper/tex/ComputGeosci/figures/QualitativelyDifferentSameGCD.pdf}
    \caption{}\label{fig:QualitativelyDifferentSameGCD}
  \end{subfigure}
  \begin{subfigure}[htbp]{1.01\textwidth}
    \centering
    \includegraphics[width=1.01\linewidth]{../../paper/tex/ComputGeosci/figures/GCDlargeButIndistinguishable.pdf}
    \caption{}\label{fig:GCDlargeButIndistinguishable}
  \end{subfigure}
  \caption[Examples showing great circle distance is a bad index of
  similarity]{(a) How the
  average GCD similarity metric ignores path geometry: \emph{Pair}\textbf{1}
  (circles and squares, left) is clearly more similar than \emph{Pair}\textbf{2}
  (circles and triangles, right), but for both pairs each GCD remains constant.
  (b) How GCD also ignores spatial uncertainties. The average GCD separation
  between coeval poles is smaller for \emph{Pair}\textbf{1} (circles and
  squares, left) than \emph{Pair}\textbf{2} (circles and triangles, right). But
  if spatial uncertainties (plotted as 95\% confidence ellipses) are considered,
  this ranking is not trustworthy: it is \emph{Pair}\textbf{2} that is
  statistically indistinguishable from the reference path. Azimuthal
  Orthographic projection.}\phantomsection\label{fig:GCDbadIndex}
\end{figure*}

We have developed an improved quantitative method of calculating the similarity
between two APWPs, or coeval segments of APWPs, in the form of a composite
difference score that compares both their spatial overlap and geometry. Our
method incorporates statistical significance testing, allowing paths with
associated spatial uncertainties to be rigorously compared to each other. The
validity and effectiveness of this method, and its superior discrimination
compared to a mean GCD score, are demonstrated by comparing the published APWP
of the North America Plate to seven derivative paths with different degrees of
spatial and geometric noise applied (Fig.~\ref{fig:2traj}).

We also test our algorithm on real paleomagnetic data, demonstrating that this
tool can be used to quantitatively assess the effects of different corrections
(in this case, bulk corrections for inclination shallowing in sediments) on the
similarity between APWPs from different continents.

\section{Methods}

\subsection{Comparing Apparent Polar Wander Paths (APWPs)}

An APWP consists of a sequence of ($\mathbf{n}$) mean poles,
$\mathbf{P_1,P_2,\ldots,P_n}$, which average the published paleopoles from a
particular continent for a particular time interval. Each mean pole has
associated longitude ($\phi$), latitude ($\lambda$), and age ($\mathbf{t}$).
Spatial uncertainty is represented by a 95\% confidence ellipse described by
semi-major axis $\mathbf{dm}$ with azimuth $\beta$ (angle east of north) and
perpendicular semi-minor axis $\mathbf{dp}$ (e.g. Fig.~\ref{fig:2traj}).

For two continents that were once part of a supercontinent, their APWPs for the
period should perfectly overlap when rotated into a common reference frame
(Fig.~\ref{fig:nambal_same_geometry}). However, due to (i) the spatial
uncertainty associated with the mean poles that form an APWP, (ii) differences
in the density and quality of data available to calculate paleopoles for
different continents in coeval time periods, and (iii) uncertainties and
possible errors in the rotations used to represent past plate motions, a perfect
match will not be obtained. Instead, two `matching' paths should share a
generally similar geometry, and largely overlap with each other when rotated
into a common reference frame. A quantitative measure of the spatial and
geometric similarity between these two paths should ideally allow us to
distinguish between non-identical paths that are similar within the associated
uncertainties, and non-identical paths that are actually different, due to
differential motion between the two continents or a poorly constrained
reconstruction.

\begin{figure*}[tbp]
  \captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj1}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairA.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj2}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairB.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj3}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairC.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj4}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairD.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj5}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairH.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj6}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairE.pdf}
  \end{subfigure}
\end{figure*}%
\begin{figure*}[!ht]\ContinuedFloat\captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:2traj7}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/pairF.pdf}
  \end{subfigure}
  \caption[Seven examples of APWP comparisons]{APWP pairs used to validate new
path comparison method. In each case the Phanerozoic APWP for Laurentia/North
America (squares, bold line) at 10 million-year (Myr) timesteps~\cite[``RM''
column of its Table 3]{T12}, is compared to a transformed copy (circles, dashed
line): (a) 1\degree\ finite rotation applied to all the mean poles and their
95\% uncertainty ellipses around an Euler pole at (125\degree\ E, 88.5\degree\
S); (b) as (a), but 15\degree\ rotation around same Euler pole; (c) after
rotation as in (b), the orientation of each APWP segment is randomised whilst
keeping their GCD length and the coeval poles' GCD fixed; (d) after rotation as
in (b), the bearing between coeval poles is randommised whilst keeping their GCD
spacing fixed; (e) as (c), but with randomly varied and relatively smaller
associated spatial uncertainty; (f) as (a), but 45\degree\ rotation around same
Euler pole; (g) as (f), but with zero associated spatial uncertainty. Azimuthal
Orthographic projection.}\phantomsection\label{fig:2traj}
\end{figure*}

\subsection{APWP Pairs Used in This Study}

To assess the performance of the evaluation method developed here, we apply it
to seven different path scenarios (Fig.~\ref{fig:2traj}) generated from
transformations of the 530\textendash0 Ma Phanerozoic APWP for
Laurentia~\citep{T12}. Almost exactly identical paths generated by rotating one
by a degree around an Euler pole at (-55\degree, 88.5\degree)
(\emph{Pair}~\textbf{\subref{fig:2traj1}}, Fig.~\ref{fig:2traj1}) represent an
ideal case of matching paths in the same spatial reference frame. Matching paths
that have been rotated out of the same reference frame by small (15\degree;
\emph{Pair}~\textbf{\subref{fig:2traj2}}, Fig.~\ref{fig:2traj2}) and large
(45\degree, \emph{Pair}~\textbf{\subref{fig:2traj6}}, Fig.~\ref{fig:2traj6})
amounts around the same rotation pole represent small and large reconstruction
errors, respectively. Random noise added to the path
(\emph{Pair}~\textbf{\subref{fig:2traj3}}, Fig.~\ref{fig:2traj3}, and
\emph{Pair}~\textbf{\subref{fig:2traj4}}, Fig.~\ref{fig:2traj4}) or the
associated uncertainties (\emph{Pair}~\textbf{\subref{fig:2traj5}},
Fig.~\ref{fig:2traj5}) represent differences in data source and/or quality. In
the final pair (\emph{Pair}~\textbf{\subref{fig:2traj7}},
Fig.~\ref{fig:2traj7}), spatial noise has been reduced by removing
\emph{Pair}~\textbf{\subref{fig:2traj6}}'s pole uncertainties.

These seven cases allow evaluation of the performance of any path comparison
metric across a range of different spatial and geometric similarities. To be
successful, such a metric must distinguish pairs with high spatial and geometric
similarity (\emph{Pair}~\textbf{\subref{fig:2traj1}}) from pairs with lower
spatial (\emph{Pair}~\textbf{\subref{fig:2traj2}},
\emph{Pair}~\textbf{\subref{fig:2traj6}},
\emph{Pair}~\textbf{\subref{fig:2traj7}})
or geometric (\emph{Pair}~\textbf{\subref{fig:2traj3}},
\emph{Pair}~\textbf{\subref{fig:2traj4}}) or both
(\emph{Pair}~\textbf{\subref{fig:2traj5}}) similarities.

To achieve more robust discrimination than the mean GCD, we propose combining a
metric for spatial misfit (Significant Spatial Difference) with metrics for
geometric difference (Significant Length Difference and Significant Angular
Difference) using a weighted linear summation, as described in the following
sections.

\subsection{Significant Spatial Difference}\label{sec:sigDs}

As in previous quantitative comparisons~\citep[for example]{S07,T08}, the
spatial separation of two APWPs is defined by the average GCD distance between
their coeval poles, but we add a filter for spatial uncertainty based on the
bootstrap approach~\citep{T91}. 1000 bootstrapped mean directions for each pole
in a coeval pair were generated (the exact sampling method is dependent on the
available information for the pole\textemdash{}see Appendix~\ref{appen4chp2} for
a full description) and their cumulative distributions in Cartesian coordinates
were compared~\citep{T91}. Pairs that could not be distinguished at the 95\%
confidence interval had their GCD separation set to 0 prior to calculation of
the mean GCD distance for all pairs. This distance is then normalised by
dividing by 50\degree, which is referred to the empirical fact that a 95\%
confidence ellipse major semi-axis of about 25\degree\ is considered
unacceptably large by paleomagnetists~\citep{B92}, to obtain the significant
spatial difference $d_s$. A $d_s$ of zero indicates that the two paths are
statistically indistinguishable from each other.

Using the same definition of $d_s$ above except testing if coeval poles are
distinguishable, the spatial difference $\Delta_s$ is defined as a comparison to
help analyze $d_s$ in the following real-data applications.

\subsection{Shape Difference}

The shape of an APWP is determined by the orientations and lengths of its
geodesic segments, which are related to the location of the Euler stage pole
that describes plate motions, and the rotation rate about that pole,
respectively. The geometric similarity of two APWPs can therefore be assessed
by comparing (i) the bearings, and (ii) the lengths of their coeval segments
(Fig.~\ref{fig:direcdiff}), with the assumption that similar geometries are
generated by a common set of stage rotations.

\begin{figure}[tbp]
  \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/dg1.pdf}
  \caption[Geometric difference definition between two APWPs]{Geometric
    differences between coeval sections of two different APWPs ({\bf
    seg$_1^I$-seg$_2^I$-seg$_3^I$} \& {\bf seg$_1^{II}$-seg$_2^{II}$-seg$_3^{II}$})
    can be described by comparing segment lengths (e.g.\ {\bf l$_2^I$} versus
    {\bf l$_2^{II}$}) or changes in bearing of coeval segments relative to their
    previous segment (e.g.\ {\bf $\alpha_3^I$} versus {\bf $\alpha_3^{II}$}).
    Segments are along great circles (blue dashed lines). Azimuthal Orthographic
    projection.}\label{fig:direcdiff}
\end{figure}

\subsubsection{Length Difference}
The length difference between the two APWPs $traj^I$ and $traj^{II}$ is the
absolute sum of differences between the lengths of coeval path segments (e.g.
$l_2^I$ versus $l_2^{II}$, $l_3^I$ versus $l_3^{II}$, $l_4^I$ versus $l_4^{II}$,
Fig.~\ref{fig:direcdiff}), normalised by dividing by the possible maximum
distance the pole could wander during the whole period $t_n-t_1$, such that:
%
\begin{equation*}
  \Delta_l = \frac{\sum\limits_{k=2}^n | l_k^I - l_k^{II} |}{D_{polar}\times(t_n-t_1)} ,
  \quad\forall k \in \{2,3,\ldots,n\},
\label{eq:ld}
\end{equation*}
%
where $|l_k^I - l_k^{II}|$ is the length difference of one pair of coeval
segments for an APWP pair ($traj^I$ and $traj^{II}$), e.g. $|l_2^I - l_2^{II}|$
for the beginning coeval segment pair. The normalising parameter $D_{polar}$ is
2.7\degree/Myr, derived from estimates of magnitude of maximum plate
velocity~\cite[up to about 30 cm/year]{S09,K14}. A $\Delta_l$ approaching 1
would result from a comparison between a virtually stationary APWP and one
associated with a rapidly moving plate.

\subsubsection{Angular Difference}

The angular difference describes the degree of consistency between the
polar-wandering directions (defined as the bearing of the older pole in a
segment with respect to the younger one) of two APWPs. In order to robustly
compare two APWPs that have not necessarily been rotated into the same reference
frame, it is more useful to define the APWP geometry relative to the path
itself, rather than an external reference frame. Therefore the bearing of a
segment is expressed as the change in geographic bearing with respect to the
previous segment ($\alpha_3$ and $\alpha_4$, Fig.~\ref{fig:direcdiff}). For
example, $\alpha_3^I$ is the result of subtracting the geographic azimuth
$\theta_{2y}^I$ from $\theta_{1o}^I$, where ``y'' stands for young end of
segment and ``o'' for old end of segment. The first segment cannot record a
relative bearing change: a path with n poles therefore consists of n-1 segments
which are described by n-2 relative angles. The defined range of bearing values
is set as -180\degree\ to 180\degree, with clockwise (east) changes in direction
defined as positive, e.g. {\bf $\alpha_3^{II}$} and {\bf $\alpha_4^{II}$}, and
anticlockwise (west) changes defined as negative, e.g. {\bf $\alpha_3^I$} and
{\bf $\alpha_4^I$}.

The angular difference $\Delta_a$ between two paths $traj^I$ and $traj^{II}$ can
then be defined as
%
\begin{equation*}
  \Delta_a = \frac{\sum\limits_{k=3}^n \Delta\alpha_k}{180\times(n-2)},
\label{eq:ad}
\end{equation*}
%
where
%
\begin{equation*}
\Delta\alpha_k =
\left\{
\begin{array}{@{}ll@{}}
| \alpha_k^I - \alpha_k^{II} |, & \text{if}\ | \alpha_k^I -
\alpha_k^{II} | \leq180; \\
360 - | \alpha_k^I - \alpha_k^{II} |, &
  \text{otherwise}.\quad\quad\quad\quad\quad\forall k \in \{3,4,\ldots,n\}.
\end{array}
\right.\label{eq:diffAziChange}
\end{equation*}

$\Delta_a$ is normalised by the maximum possible angular deviation of
180\degree. A score of 0 indicates exactly matching changes in the bearing of
coeval segments along the length of the two paths, and a score of 1 indicates
all segment bearings are antiparallel.

\subsubsection{Significance Testing of Shape Difference}\label{sec:shapeSigTest}

Due to associated spatial uncertainty, the mean poles in an APWP trace out one
possible path within a range of possible geometries (Fig.~\ref{fig:Fig5a}). If
the length and angular difference scores for one path fall within the range of
possible scores for the other, two APWPs may not in fact be significantly
different from each other. Significance testing for the shape difference scores
is performed on each coeval segment pair as follows (Fig.~\ref{fig:Fig5b}):

\begin{figure*}[!ht]
  \captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:Fig5a}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/Fig5a.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:Fig5b}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/Fig5b.pdf}
  \end{subfigure}
\caption[Testing on geometry]{Significance testing for the geometric metrics,
$d_l$ and $d_a$. (a) Illustration of how paths traj$^I$ and traj$^{II}$ can be
re-sampled within their uncertainty ellipses, with B being a possible
trajectory of traj$^I$ and C being a possible trajectory of traj$^{II}$. (b)
Upper: If the 95\% confidence interval (black vertical lines are its upper and
lower bounds) for the distribution of difference scores HIST1, generated by
comparing multiple resamplings of traj$^I$ with the original trajectory (A versus B)
does not overlap with the 95\% confidence interval (bounded by white vertical
lines) for the distribution of scores HIST2, generated by comparing resamplings
of traj$^I$ and traj$^{II}$ (B versus C), then the original difference score for
traj$^I$ and traj$^{II}$ is statistically distinguishable; Lower: If the
confidence intervals overlap, then the two paths are not
distinguishable.}\label{fig:Fig5}
\end{figure*}

\begin{itemize}
\item A bootstrapped distribution of possible geometries for each segment in a
  path can be created by resampling the two mean poles that define the original
  segment, in the same manner as described in Section~\ref{sec:sigDs} and
  Appendix~\ref{appen4chp2}.
\item A histogram of statistically indistinguishable length and/or angular
  difference scores (HIST1, Fig.~\ref{fig:Fig5b}) is created by comparing the
  resampled paths with the original for each $traj^I$ segment.
\item This distribution is then compared to the histogram of difference scores
  created by resampling the coeval segments of $traj^I$ and $traj^{II}$ (HIST2,
  Fig.~\ref{fig:Fig5b}).
\item If the two bootstrapped distributions HIST2 and HIST1 do not overlap at
  the given significance level (e.g.\ the upper and lower bounds of a 95\%
  confidence intervals, Fig.~\ref{fig:Fig5b}), then the difference score is
  interpreted to be significant. If not, then the bearings or lengths of the
  coeval segments are indistinguishable.
\end{itemize}

These tests allow a filter for spatial uncertainty to be added to the $\Delta_l$
and $\Delta_a$ metrics to define significant length difference $d_l$ and
significant angular difference $d_a$: prior to summation and normalisation, the
difference score is set to zero for the coeval segments of $traj^I$ and
$traj^{II}$ that are statistically indistinguishable.

\subsection{Composite Path Difference ($\mathcal{CPD}$)}

The three difference measures described above can be combined into a composite
path difference ($\mathcal{CPD}$) by means of a simple linear weighting rule,
%
\begin{equation*}
\mathcal{CPD} = W_s \cdot d_s + W_l \cdot d_l + W_a \cdot d_a
\label{eq:cpd}
\end{equation*}
%
for $0 < W_s,W_l,W_a < 1$, where $W_s$, $W_l$ and $W_a$ are weighting
coefficients that sum to 1. Different weighting values allow the relative
influences of spatial and geometric (length and angular) similarity to be
varied (Section~\ref{sec:wDis}).

\subsection{Fit Quality (FQ)}\label{sec:FQ}

The three metrics are all tested to be significant based on the spatial
uncertainties. However, the larger the uncertainties are, the less trustworthy
the significant difference scores are (Fig.~\ref{fig:FitQ}). Accordingly, we
bring in a concept of ``Fit Quality'', along the classification scheme of the
reversal test~\citep{M90}. For each mean pole of an APWP, it assigns a score
based on the size of the spatial uncertainty (radius: A95, or
($\mathbf{dm}+\mathbf{dp}$)/2): 1 if
it is $\leq$5\degree, 2 if 5\degree$<$r$\leq$10\degree, 3 if
10\degree$<$r$\leq$20\degree, and 4 if it is $>$20\degree. These values are
averaged for each APWP to give a ``Fit Quality'' score (from 1 to 4) for the
difference score. This is then converted into an A/B/C/D letter grade, A if the
average is $<$1.5, B if 1.5$\leq$avg$<$2.5, C if 2.5$\leq$avg$<$3.5, and D if it
is $\geq$3.5, to indicate how easy it is to generate a low difference score. In
other words, an A grade indicates that most poles are well-constrained and so it
is fairly hard to have an indistinguishable path and a low difference score. A D
grade indicates that most poles have large uncertainties so it is much easier to
have a low difference score.

In addition, a short APWP segment tends to result in overlap of its two end
poles' spatial uncertainties. For example, if an APWP is generated at intervals
of 10 Myr, the longest realistic segment-length would be about
27\degree~\cite[the maximum rate of plate movement is about 30 cm/yr]{S09,K14}.
So the uncertainty size needs be less than ${\sim}13.5$\degree\ on average to
make the segment length trustworthy. Therefore, to a certain extent, the ``Fit
Quality'' also reflects the quality of the length metric if we give each path a
grade for an APWP pair, e.g. A-A. The angular metric’s quality is related to
both coeval mean poles and successive mean poles, so it has already been
involved in the spatial and length quality. Moreover, given the fact that
usually the significant length and angular differences are much lower than
the significant spatial difference (e.g.,
Fig.~\ref{fig:sspercni},~\subref{fig:sd2ni},~\subref{fig:ld2ni},~\subref{fig:ad2ni}
and Fig.~\ref{fig:LauruG},~\subref{fig:LGsS},~\subref{fig:LGlS},~\subref{fig:LGaS}),
the ``Fit Quality'' is capable to manifest the overall quality of all the three
metrics and we should trust the difference scores if we get a B-B grade at
least.

\begin{figure}[tbp]
  \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/expect.pdf}
  \caption{Difference scores from APWPs with large uncertainties are less
trustworthy.}\label{fig:FitQ}
\end{figure}

For example, \emph{Pairs}~\textbf{\subref{fig:2traj1}},
\textbf{\subref{fig:2traj2}}, \textbf{\subref{fig:2traj3}},
\textbf{\subref{fig:2traj4}} and~\textbf{\subref{fig:2traj6}}'s fit quality
score is all 1.809\textendash1.809, so their fit quality is B\textendash{}B,
which means that the mean poles in these APWP pairs have intermediate
uncertainties on average so it is relatively hard to have a low difference
score. \emph{Pair}~\textbf{\subref{fig:2traj5}}'s fit quality is B\textendash{}A
(1.809\textendash1.085). \emph{Pair}~\textbf{\subref{fig:2traj7}}'s fit quality
is A\textendash{}A (0\textendash0).

\section{Results and Discussion}

\subsection{Discrimination of Difference Metrics}

The performance of the individual metrics were tested by generating and ranking
scores for each pair in Fig.~\ref{fig:2traj}. Scores for comparisons of the full
530 Myr paths, and sequential 100\textendash130 Myr subsections were calculated
for path pairs with (Fig.~\ref{fig:ssperc}) and without
(Fig.~\ref{fig:sspercni}) poles with zero spatial uncertainties at 350 Ma, 360
Ma, 380 Ma, 390 Ma, 450 Ma, 460 Ma and 520 Ma calculated using linear
interpolation by~\citet{T12}.

\subsubsection{$d_s$}

If ordered only in terms of spatial similarity, the desired order for the seven
APWP pairs (Fig.~\ref{fig:2traj}), from most similar (lowest $d_s$) to least
similar (highest $d_s$) is
%
\begin{equation}
\emph{Pair}~\textbf{\subref{fig:2traj1}} <
\emph{Pair}~\textbf{\subref{fig:2traj2}} \approx
\emph{Pair}~\textbf{\subref{fig:2traj3}} \approx
\emph{Pair}~\textbf{\subref{fig:2traj4}} <
\emph{Pair}~\textbf{\subref{fig:2traj5}} <
\emph{Pair}~\textbf{\subref{fig:2traj6}} <
\emph{Pair}~\textbf{\subref{fig:2traj7}}.
\label{eq:expOrds}
\end{equation}

This ordering is based largely on the mean GCD separations of each pair
(Fig.~\ref{fig:2traj}), but also takes uncertainties into account: the
relatively smaller uncertainty ellipses of
\emph{Pair}~\textbf{\subref{fig:2traj5}} and
\emph{Pair}~\textbf{\subref{fig:2traj7}} should lead to a higher $d_s$ than
\emph{Pairs}~\textbf{\subref{fig:2traj2}}-\textbf{\subref{fig:2traj4}} and
\emph{Pair}~\textbf{\subref{fig:2traj6}}, respectively. Without significance
testing, $\Delta_s$ is directly proportional to mean GCD (Fig.~\ref{fig:sd1}),
which does not result in unique $\Delta_s$ for
\emph{Pair}~\textbf{\subref{fig:2traj5}} and
\emph{Pair}~\textbf{\subref{fig:2traj7}}. Significance testing reproduces the
desired order (Fig.~\ref{fig:sd2}).

Most $d_s$ scores are also reduced with significance testing, with the largest
reductions occurring where the path separations are low and 95\% confidence
ellipses for coeval poles are more likely to overlap (e.g., $d_s$ for
\emph{Pair}~\textbf{\subref{fig:2traj1}} is reduced by 85\%, $d_s$ for
\emph{Pair}~\textbf{\subref{fig:2traj6}} is reduced by $<$2\%). $d_s$ that
approach 0 for the 0\textendash100 Ma and 100\textendash200 Ma sub-paths, which
are located close to the Euler pole used to separate the pairs and therefore
remain in close proximity even after large rotations, also illustrate this
effect. With significance testing, the 0\textendash100 and 100\textendash200 Ma
sub-paths of \emph{Pair}~\textbf{\subref{fig:2traj6}} and
\emph{Pair}~\textbf{\subref{fig:2traj7}} can still be distinguished
(Fig.~\ref{fig:sd2ni} and also Fig.~\ref{fig:sd2}), due to no spatial
uncertainty assigned to \emph{Pair}~\textbf{\subref{fig:2traj7}}. In contrast,
the older 300\textendash400 Ma and 400\textendash530 Ma sub-paths have a larger
$d_s$ than the whole path. This is because the 350 Ma, 360 Ma, 380 Ma, 390 Ma,
450 Ma, 460 Ma and 520 Ma pole coordinates are interpolated~\citep{T12}, and
thus have no assigned spatial uncertainty on any of the test paths. Without the
interpolated poles, $d_s$ is always zero for
\emph{Pair}~\textbf{\subref{fig:2traj1}} and any of its sub-paths
(Fig.~\ref{fig:sd2ni}), which is expected.

Even with significance testing, $d_s$ for
\emph{Pairs}~\textbf{\subref{fig:2traj2}}-\textbf{\subref{fig:2traj4}} are the
same (Fig.~\ref{fig:ssperc},~\subref{fig:sd1},~\subref{fig:sd2}) despite their different
geometries (Fig.~\ref{fig:2traj},~\subref{fig:2traj2},~\subref{fig:2traj3}
and~\subref{fig:2traj4}), because GCDs between coeval poles and their
uncertainties are the same. This result emphasises that a spatial difference
metric alone cannot discriminate these pairs from each other. The comparison of
\emph{Pairs}~\textbf{\subref{fig:2traj3}} and \textbf{\subref{fig:2traj5}}
indicates that well-constrained mean poles with lower uncertainties make it
relatively harder to have an indistinguishable APWP and a low difference score.

\begin{figure*}[!ht]
  \captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:sd1}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_s_dif.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:sd2}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/ds.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ld1}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_l_dif.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ld2}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/dl.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ad1}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_a_dif.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ad2}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/da.pdf}
  \end{subfigure}
\caption[Spatial, length, angular differences and significant differences of
seven exaple APWP pairs]{Spatial, length and angular differences between two
paths of the seven APWP pairs shown in Fig.~\ref{fig:2traj}. Left column:
results without significance testing imposed in the metric; Right column:
results with significance testing. Note that the spatial difference results of
\emph{Pairs}~\textbf{\subref{fig:2traj2}},~\textbf{\subref{fig:2traj3}}
and~\textbf{\subref{fig:2traj4}} are always the same for both the untested case
(a) and the tested case (b). In addition, for those segments that do not begin
from 0 Ma, their beginning segments are different from the 0\textendash100 Ma
sub-path's and the full path's. For example, for the 200\textendash300 Ma
sub-path, its beginning segment is the 200\textendash210 Ma
one.}\label{fig:ssperc}
\end{figure*}

\begin{figure*}[!ht]
  \captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:sd1ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_s_difni.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:sd2ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/dsni.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ld1ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_l_difni.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ld2ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/dlni.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ad1ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/mean_a_difni.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:ad2ni}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/dani.pdf}
  \end{subfigure}
\caption[Spatial, length, angular differences and significant differences of
seven exaple APWP pairs (without interpolations)]{Spatial, length and angular
differences between two paths of the seven APWP pairs with no interpolated poles
shown in Fig.~\ref{fig:2traj}. Left column: results without significance testing
imposed in the metric; Right column: results with significance testing. See
explanation of Fig.~\ref{fig:ssperc}.}\label{fig:sspercni}
\end{figure*}

In summary, as Fig.~\ref{fig:sd2ni} and also Fig.~\ref{fig:sd2} illustrate,
$d_s$, scaling with mean significant GCD, reproduces the expected order of
spatial similarity (Order (\ref{eq:expOrds})) for the full path. It also
compensates for the deficiency of the algorithm without statistical test
(Fig.~\ref{fig:sd1ni} and Fig.~\ref{fig:sd1}) in differentiating
\emph{Pair}~\textbf{\subref{fig:2traj6}} and
\emph{Pair}~\textbf{\subref{fig:2traj7}}. Although our algorithm also works for
APWPs with interpolations (e.g., Fig.~\ref{fig:ssperc}; see how we do
significance testing on the interpolated mean poles in
Appendix~\ref{appen4chp2}), a meaningful and valid analysis should be based on
the results with uninterpolated paleomagnetic APWPs (e.g.,
Fig.~\ref{fig:sspercni}).

\subsubsection{$d_l$}

When ordered only according to the length similarity $d_l$ the expected order is
%
\begin{equation}
0 = \emph{Pair}~\textbf{\subref{fig:2traj1}} =
\emph{Pair}~\textbf{\subref{fig:2traj2}} =
\emph{Pair}~\textbf{\subref{fig:2traj3}} =
\emph{Pair}~\textbf{\subref{fig:2traj5}} =
\emph{Pair}~\textbf{\subref{fig:2traj6}} =
\emph{Pair}~\textbf{\subref{fig:2traj7}} <
\emph{Pair}~\textbf{\subref{fig:2traj4}}
\label{eq:expOrdl}
\end{equation}

Because only the path generated for \emph{Pair}~\textbf{\subref{fig:2traj4}}
allowed the length of coeval segments to vary, it is expected that other five
pairs of APWPs have zero $d_l$ for both the full-path and the five specified
sub-paths even prior to significance testing
(Fig.~\ref{fig:sspercni},~\subref{fig:ld1ni},~\subref{fig:ld2ni}), and this
expected order is trivially reproduced. The effect of significance testing
(Fig.~\ref{fig:ld2ni}) is to substantially reduce $d_l$. Many segment length
differences do not pass the significance test because the angular uncertainties
of the poles that define individual segments are large compared to the length of
those segments.

\subsubsection{$d_a$}

If ordered only according to angular similarity $d_a$, the expected order is
%
\begin{equation}
0 = \emph{Pair}~\textbf{\subref{fig:2traj1}} =
\emph{Pair}~\textbf{\subref{fig:2traj2}} =
\emph{Pair}~\textbf{\subref{fig:2traj6}} =
\emph{Pair}~\textbf{\subref{fig:2traj7}} \leq
\emph{Pair}~\textbf{\subref{fig:2traj3}} \leq
\emph{Pair}~\textbf{\subref{fig:2traj4}} \stackrel{?}{<}
\emph{Pair}~\textbf{\subref{fig:2traj5}}
\label{eq:expOrda}
\end{equation}

Because path geometries are not altered by a simple Euler rotation, only
\emph{Pairs}~\textbf{\subref{fig:2traj3}}, \textbf{\subref{fig:2traj4}} and
\textbf{\subref{fig:2traj5}} are expected to have a non-zero $d_a$. $d_a$ for
\emph{Pair}~\textbf{\subref{fig:2traj4}} and
\emph{Pair}~\textbf{\subref{fig:2traj5}} should be larger than
\emph{Pair}~\textbf{\subref{fig:2traj3}}'s due to more geometric variation and
lower spatial uncertainties, respectively, although the expected ordering of
\emph{Pair}~\textbf{\subref{fig:2traj4}} and
\emph{Pair}~\textbf{\subref{fig:2traj5}} is less immediately obvious from visual
inspection.

Without significance testing, non-zero $\Delta_a$ for
\emph{Pairs}~\textbf{\subref{fig:2traj3}}, \textbf{\subref{fig:2traj4}} and
\textbf{\subref{fig:2traj5}} are consistently generated for both the full path
and sub-paths (Fig.~\ref{fig:ad1ni}). $\Delta_a$ for
\emph{Pair}~\textbf{\subref{fig:2traj4}} is usually higher
(Fig.~\ref{fig:ad1ni}), but there is no discrimination between
\emph{Pairs}~\textbf{\subref{fig:2traj3}} and \textbf{\subref{fig:2traj5}},
which have the same score because geometrically they are identical. When
significance testing is applied $d_a$ is markedly reduced
(Fig.~\ref{fig:sspercni},~\subref{fig:ad2ni} versus~\subref{fig:ad1ni}), and is
actually reduced to 0 for the two youngest and oldest sub-paths in all cases.
This is somewhat expected because the segment lengths of the APWPs being tested
are often of the same order as the angular uncertainty in their spatial
position. As a result, a large range of different path geometries are possible
within the specified uncertainty bounds, and the bearing of coeval segments has
to be very large for the difference to be significant.

For the full paths and the 200\textendash300 and 300\textendash400 Ma sub-paths
where $d_a$ after significance testing is non-zero for
\emph{Pairs}~\textbf{\subref{fig:2traj4}}, \textbf{\subref{fig:2traj5}} and
(usually) \textbf{\subref{fig:2traj3}}, \emph{Pair}~\textbf{\subref{fig:2traj5}}
can now be discriminated from \emph{Pair}~\textbf{\subref{fig:2traj3}}, and
consistently has the highest $d_a$ of the 3 pairs.

In summary, our angular difference algorithm with statistical test
(Fig.~\ref{fig:ad2ni}) reproduces the expected order of angular similarity
(Order (\ref{eq:expOrda})).

\subsubsection{$\mathcal{CPD}$}

When the seven different APWP pairs (Fig.~\ref{fig:2traj}) are rank-ordered in
terms of the three criteria combined, their expected order is
%
\begin{equation}
  \emph{Pair}~\textbf{\subref{fig:2traj1}} <
  \emph{Pair}~\textbf{\subref{fig:2traj2}} \leq
  \emph{Pair}~\textbf{\subref{fig:2traj3}} <
  \emph{Pair}~\textbf{\subref{fig:2traj4}} \stackrel{?}{<}
  \emph{Pair}~\textbf{\subref{fig:2traj5}} \stackrel{?}{<}
  \emph{Pair}~\textbf{\subref{fig:2traj6}} <
  \emph{Pair}~\textbf{\subref{fig:2traj7}}
\label{eq:expOrd}
\end{equation}

In order to be useful, a path difference metric needs to reproduce this order.
Note that a question mark is put on top of the ``less than'' symbols between
\emph{Pairs}~\textbf{\subref{fig:2traj4}} and~\textbf{\subref{fig:2traj5}}, and
\emph{Pairs}~\textbf{\subref{fig:2traj5}} and~\textbf{\subref{fig:2traj6}}
because when comparing pairs with different spatial separation, geometric
difference, and relative spatial uncertainty, it can be hard to objectively
define which ``should'' have the highest similarity, and the ordering will
depend on the relative weighting of $d_s$, $d_l$ and $d_a$. If the weightings
are equal (i.e, $W_s=W_l=W_a=\frac{1}{3}$), significant $\mathcal{CPD}$ scores
for paths without interpolated poles (i.e., using scores from
Fig.~\ref{fig:sspercni}) reproduce the expected order:
%
\begin{equation*}
  \textbf{\subref{fig:2traj1}} (0) < \textbf{\subref{fig:2traj2}} (0.036) <
  \textbf{\subref{fig:2traj3}} (0.043) < \textbf{\subref{fig:2traj4}} (0.067)
  \approx \textbf{\subref{fig:2traj5}} (0.067) <
  \textbf{\subref{fig:2traj6}} (0.174) < \textbf{\subref{fig:2traj7}} (0.177),
\end{equation*}
%
however \emph{Pairs}~\textbf{\subref{fig:2traj4}}
and~\textbf{\subref{fig:2traj5}} have almost identical scores and are not
discriminated. However, their fit quality (B-B for
\emph{Pair}~\textbf{\subref{fig:2traj4}} and B-A for
\emph{Pair}~\textbf{\subref{fig:2traj5}};
Fig.~\ref{fig:sspercni},~\subref{fig:sd2ni},~\subref{fig:ld2ni}
and~\subref{fig:ad2ni}) indicates that
\emph{Pair}~\textbf{\subref{fig:2traj5}}'s $\mathcal{CPD}$ is relatively more
trustworthy. This order might also not be preserved with different applied
weights. The impact of weighting will be discussed in the following section.

\subsection{A Discussion on Weights}\label{sec:wDis}

Although $W_s$, $W_l$ and $W_a$ can be defined by user, this is subjective.

However, we do explicitly know that: when comparing two APWP pairs, our aim is
to find the one whose similarity ranks higher. A simple subtraction between
unsolved (because of unknown weights) $\mathcal{CPD}s$ can help determine
which pair's similarity ranks higher. If a positive difference is obtained, no
matter what $W_s$, $W_l$ and $W_a$ values are assigned, the subtrahend pair's
similarity ranks higher; if the difference is always negative, the minuend
pair's similarity is always higher. In addition, the difference could be always
zero. Interpretation is straightforward in these three scenarios. However, for
some pairs, a positive, zero or negative $\mathcal{CPD}$ difference could result
depending on the chosen weightings.

For example, for the full (i.e., 0\textendash530 Ma) path with no interpolated
poles, the significant spatial, length and angular differences $d_s$, $d_l$ and
$d_a$ are known (Fig.~\ref{fig:sspercni},~\subref{fig:sd2ni},~\subref{fig:ld2ni}
and~\subref{fig:ad2ni}). Also we know $W_l=1-W_s-W_a$. Then we do subtractions of
$\mathcal{CPD}s$ from each two APWP pairs:
%
\setcounter{equation}{0} %reset the equation counter
\begin{equation}
\left\{
\begin{array}{lllllllllllllll}
D_{full}^{b-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.109 W_s \\
D_{full}^{c-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.109 W_s + 0.02 W_a \\
D_{full}^{d-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.038 W_s - 0.029 W_a + 0.061 \\
D_{full}^{e-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.153 W_s + 0.047 W_a \\
D_{full}^{f-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.522 W_s \\
D_{full}^{g-a} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj1}}}
= 0.532 W_s \\
D_{full}^{c-b} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}}
= 0.02 W_a \\
D_{full}^{d-b} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}}
= -0.061 W_s - 0.029 W_a + 0.061 \\
D_{full}^{e-b} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}}
= 0.044 W_s + 0.047 W_a \\
D_{full}^{f-b} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}}
= 0.413 W_s \\
D_{full}^{g-b} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj2}}}
= 0.423 W_s \\
D_{full}^{d-c} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}}
= -0.061 W_s - 0.049 W_a + 0.061 \\
D_{full}^{e-c} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}}
= 0.044 W_s + 0.027 W_a \\
D_{full}^{f-c} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}}
= 0.413 W_s - 0.02 W_a \\
D_{full}^{g-c} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj3}}}
= 0.423 W_s - 0.02 W_a \\
D_{full}^{e-d} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}}
= 0.105 W_s + 0.076 W_a - 0.061 \\
D_{full}^{f-d} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}}
= 0.474 W_s + 0.029 W_a - 0.061 \\
D_{full}^{g-d} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj4}}}
= 0.484 W_s + 0.029 W_a - 0.061 \\
D_{full}^{f-e} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}}
= 0.369 W_s - 0.047 W_a \\
D_{full}^{g-e} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj5}}}
= 0.379 W_s - 0.047 W_a \\
D_{full}^{g-f} = \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj7}}} - \mathcal{CPD}_{full}^{\textbf{\subref{fig:2traj6}}}
= 0.01 W_s,
\end{array}
\right.\label{eq:equationSet}
\end{equation}
%
and we also have the following constraints of feasible regions
%
\begin{equation}
\left\{
\begin{array}{lll}
0 < W_s < 1 \\
0 < W_a < 1 \\
0 < W_s + W_a < 1.
\end{array}
\right.\label{eq:equationSet2}
\end{equation}

The linear equations (\ref{eq:equationSet}) subject to (\ref{eq:equationSet2})
can be graphed in the three-variable ($W_s\textrm{--}W_a\textrm{--}D$)
coordinate system (Fig.~\ref{fig:pair-cmp}). For all possible combinations of
$W_s$ and $W_a$, there is a consistent ordering of $\mathcal{CPD}$ scores such
that \textbf{\subref{fig:2traj1}} $<$ \textbf{\subref{fig:2traj2}} $<$
\textbf{\subref{fig:2traj3}} $<$ \textbf{\subref{fig:2traj4}},
\textbf{\subref{fig:2traj1}} $<$ \textbf{\subref{fig:2traj2}} $<$
\textbf{\subref{fig:2traj3}} $<$ \textbf{\subref{fig:2traj5}}, and
\textbf{\subref{fig:2traj1}} $<$ \textbf{\subref{fig:2traj2}} $<$
\textbf{\subref{fig:2traj6}} $<$ \textbf{\subref{fig:2traj7}}
(Fig.~\ref{fig:pair-cmp}). However, the ranking for
\emph{Pairs}~\textbf{\subref{fig:2traj6}},~\textbf{\subref{fig:2traj7}} and
\emph{Pairs}~\textbf{\subref{fig:2traj3}},~\textbf{\subref{fig:2traj4}},~\textbf{\subref{fig:2traj5}},
or \emph{Pair}~\textbf{\subref{fig:2traj4}} and
\emph{Pair}~\textbf{\subref{fig:2traj5}} has multiple possibilities, because
their differences can be positive, negative or zero (Fig.~\ref{fig:pair-cmp}).
For this situation, assigning equal weights is recommended (giving centroid of
all possible $D$s, Fig.~\ref{fig:pair-cmp}; see also Appendix~\ref{appen4chp2}
for testing equally likely random weights) for deciding the rank order. With
equal weights used, the order from most similar pair to least similar pair is
\textbf{\subref{fig:2traj1}} $<$ \textbf{\subref{fig:2traj2}} $<$
\textbf{\subref{fig:2traj3}} $<$ \textbf{\subref{fig:2traj4}} $\approx$
\textbf{\subref{fig:2traj5}} $<$ \textbf{\subref{fig:2traj6}} $<$
\textbf{\subref{fig:2traj7}}.
These conclusions do not contradict the expected Order (\ref{eq:expOrd}).

\begin{figure}[tbp]
  \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/xyz1.pdf}
\end{figure}
\begin{figure}[!ht]
  \addtocounter{figure}{1}
  \ContinuedFloat\caption[Criteria of pair comparisons]{Graphical depiction of $\mathcal{CPD}$
    differences ($D$) between the seven APWP pairs for full-path (0\textendash530
    Ma) comparisons. If the planes derived from the equations intersect the $D=0$
    plane at a point or in a straight line, that point or the infinite number of
    points (i.e., sets of $W_s$, $W_a$ values) on the line of intersection represent
    that the similarities of the minuend pair and the subtrahend pair are equal to
    each other. If $D>0$ or $D<0$ on the planes of the equations, the subtrahend
    pair or the minuend pair respectively owns higher similarities. The square dot
    locates the result when $W_s=W_l=W_a=\frac{1}{3}$.}\label{fig:pair-cmp}
\end{figure}

In summary, as Fig.~\ref{fig:2traj} illustrates, mean GCD has trouble
discriminating between \emph{Pairs}~\textbf{\subref{fig:2traj6}}
and~\textbf{\subref{fig:2traj7}}, and \emph{Pairs}~\textbf{\subref{fig:2traj3}}
and~\textbf{\subref{fig:2traj5}}, and also between intermediate
similarities where the differences are mainly in path geometry
(\emph{Pairs}~\textbf{\subref{fig:2traj2}},~\textbf{\subref{fig:2traj3}},~\textbf{\subref{fig:2traj4}}).
Our algorithm provides an improved solution for this problem. Obtaining
similarity order can be straightforward, such as for
\emph{Pairs}~\textbf{\subref{fig:2traj1}}, \textbf{\subref{fig:2traj2}},
\textbf{\subref{fig:2traj3}}, and \textbf{\subref{fig:2traj4}},
\emph{Pairs}~\textbf{\subref{fig:2traj1}}, \textbf{\subref{fig:2traj2}},
\textbf{\subref{fig:2traj3}}, and \textbf{\subref{fig:2traj5}}, or
\emph{Pairs}~\textbf{\subref{fig:2traj1}}, \textbf{\subref{fig:2traj2}},
\textbf{\subref{fig:2traj6}}, and \textbf{\subref{fig:2traj7}}. In other words,
when one APWP's three individual metrics are all greater than or equal to, or
less than or equal to the other one's, weighting is irrelevant. However, when
the ranking of individual metrics for a pair are not consistent (e.g.,
\emph{Pair}~\textbf{\subref{fig:2traj6}} and
\emph{Pair}~\textbf{\subref{fig:2traj3}}; Fig.~\ref{fig:pair-cmp}), obtaining
similarity order is less straightforward. When this occurs, equally weighting is
recommended for concluding the final rank order.

\subsection{Application to Real Paleomagnetic Data}

To illustrate how these metrics might be useful when applied to real
paleomagnetic data, we compare 320\textendash0 Ma APWPs for Gondwana and
Laurussia calculated using a running mean method by~\citet{T12}, with paleopoles
from sedimentary rocks both uncorrected (Fig.~\ref{fig:T12Fig13a}a) and bulk
corrected for inclination shallowing (f=0.6; Fig.~\ref{fig:T12Fig13a}b).

\begin{figure}[tbp]
  \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/not_rot_yet3.pdf}
\end{figure}
\begin{figure}[!ht]
  \addtocounter{figure}{1}
  \ContinuedFloat\caption[Reproducing Torsvik et al. 2012 Earth Sci Rev
    Fig.13a]{(a) 320\textendash0 Ma APWPs (10 Myr step) for Gondwana and
    Laurussia~\cite[rotated to Southern Africa frame using the rotations
    from]{T12}; (b) as (a), but both paths corrected for inclination shallowing;
    (c) and (d) as (a), but only Laurussia path and only Gondwana path
    respectively corrected for inclination shallowing. Note that all the
    paleomagnetic APWPs are reproduced using the same moving average method and
    same paleopoles for the APWPs in Figure 13(a) of~\citet{T12}. Azimuthal
    Orthographic projection.}\label{fig:T12Fig13a}
\end{figure}

\begin{figure*}[!ht]
  \captionsetup[subfigure]{singlelinecheck=off,justification=raggedright,aboveskip=-6pt,belowskip=-6pt}
  \centering
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGs}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320s.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGsS}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320sSig.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGl}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320l.pdf}
  \end{subfigure}
  \vspace{.5em}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGlS}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320lSig.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGa}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320a.pdf}
  \end{subfigure}
  \begin{subfigure}[htbp]{.495\textwidth}
    \centering
    \caption{}\label{fig:LGaS}
    \includegraphics[width=1\linewidth]{../../paper/tex/ComputGeosci/figures/LauruGond320aSig.pdf}
  \end{subfigure}
\caption[Spatial, length and angular differences between two paths of the four
APWP pairs shown in Fig.~\ref{fig:T12Fig13a}]{Spatial, length and angular differences between two paths of the four
APWP pairs shown in Fig.~\ref{fig:T12Fig13a}. Left column: results without
significance testing imposed in the metric; Right column: results with testing.
Note that Pair~\ref{fig:T12Fig13a}b with both APWPs corrected for inclination
shallowing is not the most similar pair according to both the untested (left
column) and tested (right column) results.}\label{fig:LauruG}
\end{figure*}

When comparing the 320\textendash0 Ma paths~\cite[in their Figure 13(a)]{T12}
observed that a bulk correction for inclination shallowing applied to poles from
sedimentary rocks reduced the mean GCD separation between poles, particularly in
the Permian section of the path.

The full-path (320\textendash0 Ma) $d_s$ scores of
\emph{Pair}~\ref{fig:T12Fig13a}a and \emph{Pair}~\ref{fig:T12Fig13a}b
(Fig.~\ref{fig:LauruG}) confirm that the corrected Gondwana and Laurussia
APWPs (\emph{Pair}~\ref{fig:T12Fig13a}b) are more similar than the uncorrected
pair (\emph{Pair}~\ref{fig:T12Fig13a}a). This difference is significant
(Fig.~\ref{fig:LGsS}), and it is principally the result of an improved fit
(lower $d_s$) in the Permian (300\textendash250 Ma) and Carboniferous to
Triassic (320\textendash200 Ma) sections. When geometry is considered,
$\Delta_l$ and $\Delta_a$ scores without significance testing
(Fig.~\ref{fig:LauruG},~\subref{fig:LGl},~\subref{fig:LGa}) are actually higher
for corrected (\emph{Pair}~\ref{fig:T12Fig13a}b) than uncorrected pair
(\emph{Pair}~\ref{fig:T12Fig13a}a), particularly in the same Permo-Triassic
segment; however, none of these differences are statistically significant
(Fig.~\ref{fig:LauruG},~\subref{fig:LGlS},~\subref{fig:LGaS}). Therefore, the
equally weighted $\mathcal{CPD}$ score is actually worse for corrected pair
(\emph{Pair}~\ref{fig:T12Fig13a}b) (0.213 versus 0.1998 for
\emph{Pair}~\ref{fig:T12Fig13a}a), whereas with significance testing applied the
$\mathcal{CPD}$ of the corrected pair is a clear improvement (0.0147 versus
0.0218). This emphasises the importance of significance testing for the
geometric scores, particularly where the spatial uncertainties are relatively
large compared to the step length, as is indicated by the fit quality score
(B-A, Gondwana-Laurussia).

Furthermore, if this analysis is extended to compare an Laurussia APWP corrected
for inclination flattening with a Gondwana APWP that has not
(\emph{Pair}~\ref{fig:T12Fig13a}c), and vice versa
(\emph{Pair}~\ref{fig:T12Fig13a}d), removing the blanket correction from
Gondwana poles has only a minor effect on the overall $d_s$, whilst removing it
from the Laurussia poles actually improves $d_s$ after significance testing
(Fig.~\ref{fig:LauruG},~\subref{fig:LGs},~\subref{fig:LGsS}). Comparison of the changes in
sub-path scores for \emph{Pairs}~\ref{fig:T12Fig13a}b,~\ref{fig:T12Fig13a}c
and~\ref{fig:T12Fig13a}d suggests that the effect of the bulk flattening
correction is sometimes positive and sometimes negative for different time
periods, supporting arguments that flattening corrections need to be more
judiciously applied~\cite[for example]{B10,B16}. A more detailed study might
further constrain the regions, continents and/or time periods where a correction
is appropriate, and those where it is not. But for the purposes of this chapter,
this overview is sufficient to demonstrate the potential usefulness of our
difference metrics when considering the effect of different techniques and
corrections used to generate an APWP\@.

\section{Conclusions}
A new synthetic evaluation method is proposed in this chapter to serve as a
numerical tool for the purpose of quantitatively matching paleomagnetic APWPs.
Multidimensional information tested by bootstrapping, such as overlap of coeval
poles and shape of paths, are taken into account in the algorithm. This method
can also be utilized to detect APWP subsections' degree of similarity by
changing trajectory beginning and end poles. As an example of how this method
can be applied, we confirm a previously published suggestion that applying a
blanket correction for inclination shallowing in sedimentary rocks does
significantly improve the fit between Carboniferous to Recent APWPs for Gondwana
and Laurussia. However a more detailed analysis also indicates that such blanket
inclination corrections are unlikely to produce the best possible fit.

\section*{Computer Code Availability}
\paragraph{Name of Code} Spherical Path Comparison (spComparison)
\paragraph{E-mail} cfu3@kent.edu
\paragraph{Year First Available} 2019
\paragraph{Hardware Required} Intel (R) Core (TM) i7\textendash6700 CPU @
3.40GHz or higher; 8 GB DDR3 RAM or higher
\paragraph{Software Required} GMT5 or higher; Python 3.6 or higher; Bash 4.4.23
or higher; Linux as the best platform, macOS also fine, for Windows further
setup needed
\paragraph{Program Language} Python 3 and Shell Scripting (Bash)
\paragraph{Program Size} 75 KB (1500 lines)
\paragraph{Details on How to Access the Source Code} The source code can be
accessed from \url{https://github.com/f-i/Spherical_Path_Comparison}. Please use
the provided Jupyter Notebook file ``demo.ipynb'' to reproduce some calculations
shown in this chapter.