Testing the Reliability of Sedimentary Paleomagnetic Datasets for Paleogeographic Reconstructions

Paleogeographic reconstructions largely rely on paleomagnetic data, mostly in the form of paleomagnetic poles. Compilations of poles are used to determine so called apparent polar wander paths (APWPs), which capture the motion through time of a particular location with respect to an absolute reference frame such as the Earth’s spin axis. Paleomagnetic datasets from sedimentary rocks are particularly relevant, because of their spatial distribution and temporal continuity. Several criteria have been proposed through the years to assess the reliability of paleomagnetic datasets. Among these, the latitudinal-dependent elongation of a given paleomagnetic directions distribution, predicted by a widely accepted paleosecular variations model, has been applied so far only to investigate inclination flattening commonly observed in sedimentary rocks. We show in this work that this concept can be generalized to detect “contamination” of paleomagnetic data derived from tectonic strain, which is not always detected by field observation only. After generating different sets of simulated geomagnetic directions at different latitudes, we monitored the variations in the shape of the distributions after applying deformation tensors that replicate the effect of increasing tectonic strain. We show that, in most cases, the “deformation” of the dataset can be detected by elongation vs. inclination ratios not conforming to the values predicted by the paleosecular variations model. Recently acquired paleomagnetic directions and anisotropy of magnetic susceptibility (AMS; a parameter very sensitive to tectonic strain) data from New Caledonia verifies the results of these simulations and highlights the importance of measuring AMS when using sedimentary paleomagnetic data for paleogeographic reconstruction. We suggest to include always AMS measurement and analysis of the distribution shape to assess sedimentary paleomagnetic data used for paleogeographic reconstructions.


INTRODUCTION Paleomagnetism and Paleogeographic Reconstruction
Plate tectonics is a unifying theory that provides a solid background for understanding fundamental processes occurring on Earth. In this model, the outer shell of the Earth consists of moving lithospheric plates and their past movement can be traced using geological data . The relative and absolute motion of the tectonic plates is reconstructed using combinations of ocean floor magnetic anomalies, hot-spot tracks, and paleomagnetic data. The oldest hotspot track in the South Atlantic is ∼130 Ma, and the ocean floor magnetic anomalies allow estimating the relative finite rotation of plates back to ∼180 Ma Torsvik et al., 2008;DeMets et al., 2010;Wang et al., 2019). Paleomagnetic data are thus fundamental for paleogeographic reconstructions, particularly for pre-Jurassic times.
The basic element for all paleomagnetic-based paleogeographic reconstructions is the paleomagnetic direction. Paleomagnetic directions, measured at a given locality in rocks of any age, are expressed in the form of declination and inclination (D, I). Each direction, can be converted into a virtual geomagnetic pole (VGP), which is the point on the Earth's surface where the imaginary pole that would results in the measured D and I is located. To account for paleosecular variation of the geomagnetic field, an adequate number of VGPs should be averaged to determine a paleomagnetic pole. Subsequently, compilations of consecutive poles can be combined to establish pole paths, so-called apparent polar wander paths (APWPs). APWPs are normally plotted as a series of points "wandering" away from the geographic pole with increasing age, where paleopoles based on rocks with very young ages plot close to the geographic pole. The adjective "apparent" comes from the fact that the geographic pole does not move but APWPs rather reflect the changing orientation and distance of a plate with respect to the (fixed) geographic pole (Creer et al., 1954;Van der Voo, 1993;Besse and Courtillot, 2002;Torsvik et al., 2008). At the base of this approach lies the assumption that the geomagnetic field averaged over a few thousand years can be approximated by the one generated by a geocentric axial dipole (GAD), with the characteristic that the paleomagnetic pole, obtained by averaging the available VGPs, and the geographic poles coincide (e.g., Butler, 1992;Tauxe, 2010).

Geomagnetic Paleosecular Variation
Analyses of paleomagnetic data compilations for the last 5 Myr (McElhinny and McFadden, 1997;MM97) revealed that the recent geomagnetic field is largely dominated by the GAD component (Tauxe and Kent, 2004;Tauxe, 2005). However, a snapshot of the geomagnetic field in a particular moment at a given location on the globe would result in a VGP that is deviated from the one predicted from the GAD model up to ∼10°. This is the result of the secular variation (Fleming, 1946;Lund, 2018). Starting from the work of Constable and Parker (1988), Tauxe and Kent (2004) developed a model for the paleosecular variation (TK03.GAD) designed to fit the latitudinaldependent scatter of the VGPs observed in the paleomagnetic poles of the MM97 dataset. The TK03.GAD model results in predicted paleomagnetic direction distributions that are markedly North-South elongated, with a maximum elongation at the equator gradually diminishing toward the poles (Tauxe and Kent, 2004;Tauxe, 2005). Following this model, a collection of paleomagnetic directions from a locality, spanning a time interval long enough to average out secular variations, is expected to possess a distribution with a certain degree of elongation developed parallel to the mean declination.

Paleomagnetic Earth Filters
Sedimentary rocks constitutes about 66% of the rocks exposed on the Earth's surface (Blatt and Jones, 1975). Therefore, sedimentary paleomagnetic datasets are of extreme value for paleogeographic reconstructions. However, the record of the Earth's paleomagnetic field in sedimentary rocks is often "distorted" by different natural processes including tectonic deformation, remagnetization, or compaction-induced inclination flattening. These phenomena have been defined as "Earth filters" (Tauxe, 2005). In particular, sedimentary inclination shallowing of paleomagnetic directions is a long known problem, and it is described by the tangent function introduced by King (1955): where I o is the observed remanence inclination, I f is the inclination of the inducing field, and f is the flattening factor ranging from 0 (completely flattened directions) to 1 (absence of flattening). Directional sets affected by significant inclination shallowing typically possess an elongation developed perpendicular to the average declination (Tauxe and Kent, 1984). Eq. 1 is at the base of the E/I (Elongation/ Inclination) method for correcting shallow biased paleomagnetic directions (Tauxe and Kent, 2004). A given set of flattened paleomagnetic directions is corrected by applying gradually decreasing f values to each direction. The "unflattened" average inclination is reached when the elongation of the whole distribution matches the value predicted by the TK03.GAD model. This statistical technique has been successfully applied to a wide range of sedimentary rocks (Kent and Tauxe, 2005;Krijgsman and Tauxe, 2006;Tauxe et al., 2008;Dallanave et al., 2009, Dallanave et al., 2012, Dallanave et al., 2015, Dallanave et al., 2018aKirscher et al., 2014). A reliable APWP of Adria (the African promontory) since the Permian was recently compiled systematically applying the E/I correction to new and published datasets (Muttoni et al., 2013;Muttoni and Kent, 2019). Paleomagnetic inclination flattening is not the only "Earth filter". Many sedimentary rocks are affected by tectonic strain, which is not always detectable by simple field observations, and that can deflect the characteristic remanent magnetization (ChRM) directions (Cogné and Perroud, 1985;Lowrie et al., 1986;Jackson et al., 1993;Borradaile, 1997). The finite strain state of a sedimentary rock can be qualitatively evaluated by using the anisotropy of magnetic susceptibility (AMS). By monitoring the AMS of Eocene mudstones of the Southern Pyrenean Foreland Basin, Parés et al. (1999) found that during incipient deformation, the fabric evolves through a series of stages that can be essentially synthetized in three types ( Figure 1). In the absence of deformation, only compaction acts on the sediments resulting in an oblate AMS ellipsoid with a vertical minor axis ( Figure 1A). As the deformation evolves to the weak cleavage state, the fabric changes toward a prolate shape with the major axis perpendicular to the shortening direction ( Figure 1B). At the strong cleavage state, the fabric will assume a triaxial k 1 > k 2 > k 3 form (where k 1 , k 2 , and k 3 are respectively the major, intermediate, and minor axes of the AMS tensor) with the minimum axis parallel and the maximum perpendicular to the shortening direction ( Figure 1C).
In some pioneering work, ChRM directions in sedimentary rocks have been restored with an "unstrain" strategy (Cogné andPerroud, 1985, Cogné andPerround, 1987;Cogné, 1987). In this approach the paleomagnetic directions are considered as behaving like a passive line within the rock matrix, even though in some cases, due to the nature of the rock forming particles, this approach may be an oversimplification (Kligfield et al., 1983;Borradaile, 1997). The quality and reliability of the "unstrained" dataset was assessed using basic Fisher (1953) statistical criteria, whereby increasing of the precision parameter k reflects improved clustering of the directions. It is now clear that the remanence magnetization of sedimentary rocks reflects a more complex behavior of the geomagnetic field, and using just the standard Fisher (1953) statistics to evaluate the reliability of paleomagnetic datasets is an oversimplification. Monitoring the E/I of a paleomagnetic dataset in order to detect inclination flattening has become already a standard approach. Analogously, other earth filters can lead to E/I couples of values that are "unrealistic" (i.e., departing from the values expected from the TK03.GAD field model). In this paper we show that the shape of given ChRM directions set can be used to check if finite strain has undermined the reliability of the data. Using simulated paleomagnetic data, we show that strain related deflection of ChRM directions promptly modify the distribution elongation expected from the TK03.GAD model at a given latitude. We apply this concept to assess the reliability of a recently published dataset from Eocene sedimentary rocks exposed in New Caledonia (Dallanave et al., 2020).

Data Generation
To clearly visualize the effect of inclination flattening and finite strain simulation on a theoretical distribution of directions, we first generated a purely circular distribution of 30°radius around a mean direction of D 0°and I 50°( Figure 2A). A circular distribution can be considered as a simplification of a Fisher (1953) distribution, which is defined as 1) uniformly distributed around the mean direction and 2) with a occurring frequency of directions decaying exponentially with the distance form the mean (see also Fisher et al., 1987).
To evaluate the effect of finite strain on more "realistic" directions, which reflect the secular variations of the geomagnetic field, we generated six sets of paleomagnetic directions using the tk03.py Python script compiled by Tauxe et al. (2016). The sets consist of declination-inclination-intensity triplets drawn accordingly to the TK03.GAD geomagnetic field model at specified latitudes. The six sets were generated applying latitudes of 10, 15, 20, 30, 40, and 50°N, imposing a normal polarity field. For all sets we selected a number (N) of 120 directions (Supplementary Table S1). Tauxe et al. (2008) indicated that an adequate number of directions to evaluate reliably the distribution shape is > 100. Since it is uncommon to find published direction sets that exceed this number by far, we consider N 120 realistic. We excluded all directions with geomagnetic intensity lower than 10μT from the generated datasets. This is based on the assumption that the highest harmonics of the field (i.e., with the lowest intensity) are often not recorded by sediments, and if they are, they result in transitional directions that are normally excluded by directional cutoff (analogously to what done by McElhinny and McFadden (1997), for the MM97 compilation). The obtained distributions possess k values ranging from 20 to 35 ( Figure 3; Table 1), which is in agreement with published high-quality sedimentary datasets (see e.g., compilation of Muttoni et al., 2013).

Strain Simulation and Elongation Monitoring
Paleomagnetic directions are commonly defined by D and I with unit length, so the intensity of the TK03.GAD directions was not considered further after filtering. For each generated set we In a strong cleavage state the fabric is triaxial, with the minor and major axes respectively parallel and perpendicular to the shortening direction. The effect of the fabric evolution on the paleomagnetic inclintaion is shown on top of the figure; figure drawn adapting concepts from Parés et al., 1999 andBorradaile, 1997. Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 592277 calculated the elongation E, which is defined as the v 2 /v 3 ratio, where v 2 and v 3 are the intermediate and minimum eigenvalues of the directions distribution matrix (Scheidegger, 1965; Table 1). The directions of each set are first converted to Cartesian coordinates x i with i 1, 2, 3. We simulated the effect of the deformation by applying the strain matrix S, which is a symmetric matrix of elements s i,j (where i, j 1, 2, 3, and s i,j s j,i ). The new coordinates d i of each "strained" direction are obtained by: This can also be used to simulate the effect of inclination shallowing by applying a diagonal matrix S of the form s 11 s 22 > s 33 (flattening matrix F). Note that under this condition this operation is equivalent to apply Eq. 1 with f s 33 /s 11 . FIGURE 2 | Effect of strain simulation on a purely circular distribution. (A) 1 original distribution; 2 simulated inclination flattening; 3 effect of a pure compressional shear, with shortening indicated by the red arrow; 4-6 effect of combination of pre-compressional strata tilting followed by shortening; while only northward tilting is shown in equal area projections 4 and 5, projection 6 shows the effect, after correcting for the bedding tilt, of both northward (blue) and southward (red) tilting on the same distribution; 7-9 same as 4-6 but with an initial tilting of 30°and simulated shortening directed NE-SW; F flattening simulation; S shortening; T10°and T30° 10°and 30°tilting; uT correction for strained tilt. In all insets: D declination (°); I inclination (°); v 2 and v 3 medium and minimum axes of the distribution matrix (Scheidegger, 1965); E elongation (v 2 /v 3 ). (B) Diagram of elongation (E) vs. inclination (I) of the original circular distribution (1), after inclination flattening (2), and after all the listed strain simulations; 6 and 9 are the E/I couples for the blue data plotted in the corresponding projections of panel A, while 6* and 9* refer to the red data; the shape, size, and orientation of the ellipses reflects the variations of the directions distributions as determined by the distribution matrix, while the arrows behind the ellipses are the mean declination of the directions; SDZ and FDZ shallowed and flawed distribution zone, respectively.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 592277 We monitored the effect of the strain simulation on the directions distributions following the steps described as follows.

Circular Distribution
The workflow is shown in Figure 2. We first apply an inclination flattening with f ∼0.6 (detailed values of s i,j are listed in Table 2) to the circular distribution, which can be considered a reasonable assumption for sedimentary rocks (Figures 2A, diagram 1 and 2). In nature, the flattening factor is case sensitive and can easily range from 0.4 to 0.9 (see e.g. compilation of Muttoni et al., 2013). We simulated at first the simplest scenario, by applying shortening parallel to the average declination ( Table 2) as indicated by the red arrow ( Figures 2A, diagram 3). This situation, however, is not entirely realistic, because it implies no tilting of the strata, and in nature compressional strain is typically associated with folding. For this reason, in the next step we followed the case illustrated by Borradailee (1997). At the very initial stage of the deformation, low amplitude buckles, which result in tilted strata but not pervasive deformation, are subsequently passively strained, and the strain is affecting both the paleomagnetic directions and the strata attitude. The initial tilting of the strata is unknown, so we applied two different example scenarios, with initial tilting of 10°N (and 10°S), followed by N-S shortening ( Figure 2A, diagrams 4-6), and initial tilting of 30°N (and 30°S) with shortening oriented NE-SW ( Figure 2A, diagrams 7-9). We plot the strained directions, with associated  The sets are defined by their (approximated) latitude, as illustrated in Figure 3. Dec., Inc. average declination and inclination; k, α 95 precision parameter and 95% confidence angle (Fisher, 1953); E elongation of the distribution, defined as the ratio of the intermediate and minimum eigenvalues of the directions distribution matrix (Scheidegger, 1965). All simulated sets are listed in supporting Supplementary Table S1 available online. S ij elements of the matrix of Eq. 2; F flattening matrix; C and C (45°E) shortening matrices applied to the pure circular distributions, directed respectively N-S and NE-SW; S w and S t matrices of Eqs 3, 4 simulating respectively the weak cleavage and the strong cleavage state; S w (45°E) and S t (45°E) same as before but simulating shortening oriented NE-SW.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 592277 means, in tilt corrected coordinates, and the variations of E/I ratios are monitored in Figures 2B.

TK03.GAD Distributions
Similar to what has been done for the circular distribution, we applied different straining scenarios, while monitoring the variations of E/I ratios in all cases. At first, to explore the effect on the original TK03.GAD distribution shape of pure compressional strain, we applied a series of oblate deformation tensors in the form (expressed as eigenvalues) of s 1 s 2 > s 3, with s 1 vertical. The value of the (normalized) eigenvalue of s 3 ranges from 0.9 to 0.2 ( Table 3). We then varied the orientation of the eigenvector s 3 with respect to the average declination of the distributions, ranging from 0°to 20°. Geometrically, this replicates different tectonic stresses oriented parallel to the strata but with a different azimuth (Figures 4,5). After each simulated deformation we calculate the average inclination and the elongation of the distribution for comparison with the E/I value expected from the TK03.GAD model. A more complex but more realistic scenario has been applied to the direction sets simulated for low (10°N) and mid (40°N) latitudes, with the aim of recreating the fabrics shown in Figure 1. The directions have been first flattened with f ∼0.6 ( Figures 1A). The weak cleavage and the strong cleavage fabrics are assumed to be the result of the sum of the flattening fabric and a successive strain. As the weak cleavage state is characterized by a prolate fabric with k1 > k2 k3 (Figures 1B), while the strong cleavage by a triaxial fabric with k1 > k2 > k3 ( Figures 1C), we applied to the flattened directions a deformation matrix calculated by either: or where S w and S t are the deformation matrix that simulate respectively the weak and the strong cleavage final stage, P w (weak cleavage) and P t (strong cleavage) are the (case sensitive) fabric that characterize the rocks, and F −1 is the inverse of the diagonal matrix F used to simulate the inclination shallowing ( Table 2).

Circular Distribution
The effect of inclination shallowing followed by simulated shortening is shown in Figure 2. As expected, the solely flattening results in a variation of the elongation from circular to E-W elongated ( Figures 2B). This condition can be easily restored by the site-level E/I correction (Tauxe and Kent, 2004), which gradually increase the directions inclination using Eq. 1 until they are Fisher (1953) distributed (i.e., circular around the mean direction). The simple scenario of N-S shortening without rotation (i.e., bedding tilt) results in an increase of both elongation and inclination (Figures 2A,B, distribution 3). It is interesting to observe the influence of pre-existing tilting on the parameter variation, especially the inclination. The equal area projection 6 of Figures 2A shows the result of a N-S shortening simulation of distribution initially tilted by 10°N (blue ellipsis) and 10°S (red ellipsis). The final inclination is respectively ∼15°s hallower and ∼14°steeper than the one obtained by the pure shear shortening ( Figures 2B). This is also in agreement with the behavior in natural sedimentary rocks predicted by Borradaile (1997) and schematized in Figure 1. If tilting and paleomagnetic inclination are plunging toward the shortening direction, the final inclination with respect to the bedding will be flattened, while if the bedding and directions are antipodal, the final inclination will be steepened. Changing the orientation of the S matrix, will also affect the elongation direction, as shown in diagrams 7-9 of Figure 2. The inclination is less affected than during shortening parallel to the paleomagnetic declination, but the final elongation is not perpendicular to the declination, as expected if the inclination shallowing was the only filter acting on the data (Figures 2B, ellipses 9 and 9*). This concept is important also to help identifying the possible effect of finite strain in natural sedimentary rocks, as discussed below.

TK03.GAD Distributions
Effect of Strain Without Tilting Figure 3 shows the simulated undeformed distributions. The calculated E/I lies in proximity to the expected values from the TK03.GAD model. We first apply the simplified model of strain that involves no pre-deformation rigid body rotation, with initial shortening parallel to the paleomagnetic declination ( Figure 4, Group 1). The bold number within each equal area projection of Figure 4 corresponds to the applied tensor listed in Table 3. The effect on the expected E/I couplets is shown in Figure 5. The most extreme applied simulated deformations are probably very unlikely to be found in nature, especially without the folding and tilting affecting the rocks, and possibly other remagnetization effects that would modify or even replace the original directions (i.e., piezoremanent magnetization; PRM; Till et al., 2010). As expected, when s 3 is parallel with the average declination there is a progressive increase of the inclination. The minimum of the curves in Figure 5 coincides with the point where the elongation direction changes from parallel to perpendicular with respect to the declination. Before that point they are systematically below the reference TK03.GAD values. Soon after the elongation direction shifts from parallel to perpendicular with respect to the declination, the elongation increases to values higher than the reference TK03.GAD curve ( Figure 5). When increasing the angular distance between the declination and the eigenvector s 3, there is progressively less effect on the inclination. This occurs along with a stronger effect on the elongation ( Figure 5). An angle of 20°between s 3 and the distribution declination is enough to drive the E/I to unrealistic values at low levels of strain. Direction distributions at high latitude, which are predicted to be quasi-circular by the TK03.GAD model, acquire an unrealistic elongation at very low strain (Figures 5F).

Effect of Strain With Tilting
We followed a chain of events closer to reality by combining the effects on the 10°N and 40°N TK03.GAD distributions of an initial compaction flattening, followed by finite strain simulations that replicate the weak cleavage and the strong cleavage fabric observed in natural rocks. Both strains are applied to directions belonging to two limbs of a fold ( Figures 6A). We first assumed a shortening parallel to the paleomagnetic declination ( Figures 6B,C). After applying the compaction flattening matrix F ( Table 2), both the 10°N and the 40°N distributions acquire a nearly circular shape ( Figures 6C). Because the initial bedding tilt before deformation is unknown, in the case of the weak cleavage we applied an initial strata inclination of 10°, plunging both to the north (front limb) and to the south (back limb). The effect of the S w matrix of Eq. 3 on the 10°N directions set is not significant in terms of inclination, but the distributions of both limbs acquire an elongation oriented E-W and falling below the reference TK03.GAD line ( Figures 6B,C, FIGURE 4 | Effect of simulated deformations, without strata tilting, on the TK03.GAD directions datasets. The latitude of the generated dataset is also indicated, and the colors of the data are as in Figure 3. Group 1 and 2 correspond respectively to a N-S and a 20°E oriented shortening, as represented by the eigenvectors of the strain matrix (S) shown in the insets. The bold numbers inside the projections indicate the S matrix applied to the set, listed in Table 3.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 592277 distr. 1 and 2). The same strain simulation appears to have a more severe impact on the 40°N set, because both fold limbs are characterized by distribution shapes jumping above the TK03. GAD line (Fig. 6B,C, distr. 3 and 4). When applying the S t matrix (Eq. 4), which replicates the strong cleavage scenario, we assumed an initial bedding tilt of 30°for both limbs. The combined effects of stronger strain and higher initial bedding are immediately visible on the E/I ratios of both the 10°N and the 40°N sets. The inclinations of the 10°N set do not show a dramatic variation, but in the front limb case the strained distribution (which gets sensibly smaller) is strongly E-W elongated, while the back limb is still well below the TK03.GAD line ( Figures 6B,C, distr. 5 and 6).
As for the weak cleavage simulation, the 40°N distribution appears to be more affected by the S t matrix, with the front and the back limbs distributions showing respectively a marked shallowing and steepening of the inclination, parallel to a significant variation in elongation ( Figures 6B,C, distr. 7 and 8).
Applying a shortening oriented 45°clockwise, has in general a more severe impact on all distributions ( Figures 6D,E). The effect is particularly evident on the elongations; in all cases the strong cleavage scenario causes the E/I couples to assume unrealistic values. The weak cleavage simulations generate E/I values that are closer to the expected TK03.GAD reference line, nevertheless the direction of the elongation is significantly deviated (Figures 6D,E).

Strain and Distributions Shape
In the case of the circular distributions, the effect of either inclination shallowing or finite strain (or the combination of both) is readily visible, because the distribution departs from the original circular shape (E 1). With the assumed circle (30°of radius) all the simulations, except for the pure compaction shallowing, fall above what we name the "shallowed distribution zone" (SDZ). This is the area where we would find all E-W elongated ellipses affected only by compaction shallowing, and they could be restored through the "site mean" E/I correction (Tauxe and Kent, 2004). All other E/I values of this diagram (for a circle of initial 30°radius) indicate distributions that are flawed and affected by reorientation mechanisms like strain, which make them unreliable. We call this area of the diagram the "flawed distribution zone" (FDZ) ( Figures 2B). From the analysis of the circle distribution it is also evident that the angle between the maximum axis of the elongation and the mean direction declination can be used as a reliability index. In fact, inclination flattening only produces elongations that are perpendicular with respect to the declination, and thus an angle that is not perpendicular is indicating a flawed distribution. Similar concepts can be applied to the analyses of the TK03.GAD distributions. Applying the S matrix on the original TK03.GAD distributions, without imposing a prestrain rotation (i.e., tilting), has an effect that depends mostly on the angle between the shortening direction (s 3 eigenvector of the strain matrix) and the average direction declination. When they are parallel (Figure 4, Group 1) many of the E/I values of the deformed datasets, fall below the expected TK03.GAD function line, in the SDZ section of the diagram ( Figures 5A-F). Analogously to what has been described for the circular distributions, this is the area of the diagram where normally the E/I pair of sets that are affected by sedimentary inclination flattening are found. Applying the E/I unflattening method of Tauxe and Kent (2004) to such a distribution would result in a final inclination even further steepened. As the degree of the simulated strain increases, the distributions acquire unrealistic E/I pairs, falling in the FDZ. Since high latitude simulations are originally quasi circular, they are promptly deformed to unrealistic E/I values by the simulated strain ( Figures 5F). A similar effect is obtained by increasing the angle between the shortening direction and the declination (Figure 4, Group 2, Figure 5). This is the case because it affects more the distribution shape rather than the inclination. A shortening direction oriented E-W would in fact act only on the shape, which would peak straight to unrealistic values.
The effect of pre-strain tilting is affecting the final directions inclination (calculated with respect to the bedding plane) in a more complex way. With the applied boundary conditions (flattening and initial bedding tilt), especially in the strong cleavage scenario, the direction belonging to the front limb of the modeled fold tends to be flattened, while the ones belonging to the back limbs are steepened. This is in agreement with the behavior predicted by Borradaile (1997). Strain simulations with a shortening azimuth deviated from the mean declination affect more the elongation, also in terms of orientation. Any paleomagnetic set that is elongated neither parallel (as predicted by the TK03.GAD model) nor perpendicular (i.e., affected by compaction shallowing) with respect to the average declination should thus be considered potentially flawed. These analyses also show how sensitive the initial conditions like flattening and bedding tilt, which are a priori unknown, influence the final result. This is why paleomagnetic analyses for paleogeographic reconstructions should always be carried out in combination with AMS analyses, as a proxy of the presence of finite strain. Only when the AMS fabric is sedimentary, a standard E/I correction can be safely applied.

THE NEW CALEDONIA CASE
Paleomagnetic Directions and Anisotropy of Magnetic Susceptibility Dallanave et al. (2020) recently published paleomagnetic data from carbonate rocks exposed in the Koumac region of northern New Caledonia. New Caledonia is the emergent part of the northernmost Norfolk Ridge, which is in turn part of a large continental mass that is submerged for more than 90%, referred to as Zealandia (Mortimer et al., 2017;Sutherland et al., 2019). The sampled section (Sommet-Khian, 260 stratigraphic meters) consists of a massive basal pelagic micrite (0-83 m) overlain upsection by terrigenous-rich calciturbidites (83-260 m). For this work we exclude the dataset from the terrigenous-rich calciturbidite because of the general lower quality of the ChRM directions and a negative reversal test (see Dallanave et al., 2020 for details). The dataset from the massive micrite consists of 88 ChRM directions that allow a correlation of the section with Chrons C23n.2n to C21n (51-46 Ma, corresponding to the earlymiddle Eocene; Ogg, 2012). Of these 88 directions, 77 have been isolated using the interpolation proposed by Kirschvink (1980) and, to maximize the reliability, each direction was anchored (or not) to the origin of the demagnetization axes following the Bayesian criteria proposed by Heslop and Roberts (2016). Eight directions were determined by means of fisher mean on the vector end points, while three through great circles analyses (McFadden and McElhinny, 1988); (Supplementary Table S2). The dataset is statistically antipodal, with a positive bootstrap-based reversal test (Tauxe et al., 1991). Despite the high quality of the dataset, there is a mismatch between the average inclination from Sommet-Khian and the reference directions expected for this time (Veevers and Li, 1991;Torsvik et al., 2012). The average direction from New Caledonia is steeper (i.e., higher inclination) compared to the reference directions from the literature. When compared with the inclination derived from the synthetic global APWP of Torsvik et al. (2012), the resulted steepening is ∼10°. However, this reference compilation consists of a limited number of entries that has to go through multiple rotations to be plotted into Zealandia coordinates, and there is only one entry from the Australian plate (53 Ma; Torsvik et al., 2012) Table S3). The steepening decreases to 5.4°i f the inclination is compared with the APWP of Veevers and Li (1991), which is based on entries from Australia and India, but using old datasets with poorly defined quality criteria and age control ( Figures 8A).
AMS data of 89 samples from Sommet-Khian yield an oblate tensor with a low degree of anisotropy (Jelínek, 1981) (Pj 1.05). The foliation plane (defined by k 1 and k 2 axes of the AMS) does not coincide with the bedding, indicating a tectonic-induced magnetic fabric in contrast to a sedimentary fabric typically observed in absence of finite strain ( Figures 8B). The orientation of this fabric is in agreement with the general tectonic style of New Caledonia, which is the result of a Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 592277 convergence broadly oriented NE-SW, perpendicular to the length of the island (e.g., Lillie and Brothers, 1969;Cluzel et al., 2012). Tectonic strain largely influences the AMS fabric (e.g., Parés et al., 1999) but using the absolute value of the AMS to quantify the internal strain of the rock is difficult because of the number of variables that can influence it (Evans et al., 2003).  Table S3); the 90-10 Ma global synthetic apparent polar wander path from the same compilation, plotted in Australian coordinates, is also shown as thick black line; the red line with the shaded area shows that, even rotating the S.K. pole to toward reference 48 Ma one, they remain statistically distinguishable. (C) Paleogeographic reconstruction of the Southwest Pacific area in the middle Eocene (∼45 Ma); continental areas are filled in light brown except for northern Zealandia (light green); the light gray line is the Tasman Ocean spreading ridge, which activity ended during Chron 24n (∼53 Ma), and the dark gray lines are active spreading ridges; dotted black lines are the approximate position of the proto-Tonga-Kermadec subduction zone and the north-east dipping early Eocene east-dipping subduction zone inferred by the geology of New Caledonia (see text); the white star is the location of Sommet-Khian. While the reconstruction is based on the published reference data, the red line with the red shaded area is the paleolatitude and associated confidence boundaries obtained from the paleomagnetic directions from Sommet-Khian.

Analysis of the Directions Distribution From Sommet-Khian
The first step is comparing the E/I ratio of the dataset from Sommet-Khian with the one predicted by the TK03.GAD model. The original ChRM directions have an elongation that falls above the value predicted by the TK03.GAD model. Repeated calculations on 1,000 bootstrapped pseudo-distributions (Tauxe et al., 1991) estimate a (fairly large) 90% confidence boundary for the elongation that overlaps with the expected value ( Figures 8A), implying that the E/I ratio itself is, at a 90% confidence, realistic. The wide extension of the 90% confidence margins is likely due to the limited number (<100) of directions. This outlines the importance of having an adequate number of paleomagnetic directions for a fully reliable analysis. We repeated the same calculation using the dataset filtered from the three directions determined by great circle analyses and also  eliminating five directions which VGP is > 45°apart from the mean paleomagnetic pole, similarly to the cutoff applied to the MM97 dataset of McElhinny and McFadden (1997). The obtained inclination is very similar, with a slightly higher elongation ( Figures 8A). Despite the confidence margins, the fact that both datasets have an elongation higher than expected is an indication of potential contamination of the dataset, possibly steepened by the effect of strain. We tried to restore the original inclination by progressively "unstrain" the paleomagnetic directions until they match the elongation predicted by the TK03.GAD model. In order to do this we considered the bedding also as a physical entity affected by the total strain. We therefore apply the "unstrain" matrix in the form of s 1 > s 2 s 3 with progressively higher degrees of elongation, with s 1 parallel to the k 3 AMS axis (Figure 8). We applied the "unstrain" matrix to all the ChRM directions in situ (IS) coordinates as well as to the bedding plane using the eigenvalues listed in Table 4, and we monitor the E/I couplets after correcting the dataset for ("unstrained") bedding tilt ( Table 5). We applied the same process to both the whole and the filtered dataset, so that we can evaluate the effect of outlier directions within the (not-filtered) dataset. In both cases the progressive "unstrain" led to a minor increase of the inclination, which however always falls within the confidence bounds of the original distribution ( Figures 8A). The value of E changes in different fashion in the two cases: for the whole dataset it progressively decreases, but without crossing the model line. Despite that, the bootstrap-based 90% confidence interval indicates that the elongation is statistically acceptable throughout the whole "unstrain" process, and it reaches a minimum after applying matrix 4 of Table 4 ( Figures 8A).
The filtered dataset on the contrary is characterized by a steep rise of the elongation, which becomes statistically unacceptable (at a 90% confidence) after a very mild unstrain (S matrix 2 of Table 4; Figures 8A).
A close inspection of the equal area projections, plotted together with the elongation direction (Figures 8C) reveals that the elongation assumes an orientation parallel to the declination after applying unstrain matrix 6 (in the case of the whole dataset) and 1 (in the case of the filtered dataset). Having the elongation parallel to the mean declination is, together with the latitude dependent specific value of E, a condition predicted by the TK03.GAD model. The inclinations obtained are 64.1 ± 3.5°and 63.5 ± 2.9°respectively for the whole and the filtered dataset ( Figures 8A). In both cases the elongation is still statistically acceptable at a 90% of confidence. In synthesis, both datasets give similar corrected results, even if the limited number of directions results in wide elongation confidence bounds.

Paleogeographic Implications
The southwest Pacific area went through a complex tectonic and paleogeographic evolution during the Eocene, related to the inception of the Tonga-Kermadec subduction system (recent reviews can be found in Cluzel et al., 2012;Collot et al., 2020;Dallanave et al., 2018b, Dallanave et al., 2020Maurizot et al., 2020a, Maurizot et al., 2020bSutherland et al., 2020). The Tonga-Kermadec subduction initiated in the middle Eocene immediately to the east of the Norfolk Ridge, the northernmost part of which is The eigenvector are oriented as the average eigenvector of the anisotropy of magnetic susceptibility (AMS), shown in Figure 8B. The eigenvalues used for the "unstraining" process are indicated as v 1 , v 2 , and v 3 , with the orientation (V i _dec, V i _inc. declination and inclination) of the associated eigenvector. All_0 to 6 refer to the completed directional dataset (88 direction) without filtering, while Filtr_0 to 6 refer to the filtered dataset (80 directions) as specified in the text; Dec, Inc. declination and inclination; k and α 95 precision parameter and 95% of Fisher, 1953. The bedding, expressed as declination and inclination of the plunge, is also "unstrained" before applying the bedding correction to the paleomagnetic directions. The derived tilt corrected directions are given with the associated elongation (E) of the distribution and the 90% confidence boundaries determined by repeating the analysis on 1,000 bootstrapped sets of directions; (-) and (+) are respectively the lower and the higher confidence margins. The complete sets of directions are given in supporting Supplementary The history of the eastern margin of northern Zealandia before the Tonga-Kermadec subduction is still debated, because the geology of New Caledonia suggests the presence of a pre-Tonga east-dipping subduction zone (Cluzel et al., 2012;Collot et al., 2020;Maurizot, 2011;Maurizot et al., 2020a, Maurizot et al., 2020b, which is, however, difficult to reconcile with the geologic history derived from seismic profiles acquired in the Tasman area (Sutherland et al., 2017). Zealandia is separated from continental Australia by the oceanic crust of the Tasman Ocean basin, which ended its spreading activity during magnetic anomaly 24 in the early Eocene (∼53 Ma; Gaina et al., 1998;Figures 7C). Therefore, middle Eocene paleomagnetic data from Australia should be applicable also to northern Zealandia. The most widely used APWP was compiled by Torsvik et al. (2012), but a closer look at the dataset reveals that there is only one early Eocene (53 Ma) entry from Australia for the entire Paleogene based on a study on the Barrington Volcano in New South Wales (Wellman et al., 1969, later confirmed by;Hill et al., 2002). Recent authochtonous data have been published only for the Oligocene (Hansma and Tohver, 2019). The reference Australian APWP for the Paleogene (and the entire time back to the Permian) is constructed mainly using data rotated from other plates. For this reason, the paleomagnetic record of Sommet-Khian is extremely valuable, because it is the first early-middle Eocene autochthonous directional dataset from northern Zealandia (and Australia). The Eocene sedimentary record of the Koumac area lies on top of the Mesozoic basement, but New Caledonia as a whole is constituted by sedimentary, metamorphic, and igneous terranes that are overlain by allochthons of sedimentary, mafic, and ultramafic nappes obducted from the late Eocene with the opening of the North Loyalty Basin (Cluzel et al., 2001, Cluzel et al., 2005, Cluzel et al., 2006, Cluzel et al., 2012Maurizot, 2011;Maurizot et al., 2020a). The paleomagnetic declination obtained by Dallanave et al. (2020) is in general agreement with the tectonic constraints, indicating a postmiddle Eocene rotation of New Caledonia of about 60°CCW. The ChRM directions inclination is primary and should be usable (at a first inspection) as paleogeographic constraint.
The "unstrained" mean paleomagnetic inclination is ∼64°( Figures 8A). Using the dataset with the widest uncertainties (i.e., the unfiltered directions; 64.1 ± 3.5°) we obtained an average paleolatitude of deposition of 45.8°S (confidence boundaries: 50.5°S, 41.6°S; Figures 7C). We estimated precisely the expected paleolatitude of Sommet-Khian by using 14 paleomagnetic poles falling within the 53 and 43 Ma age window from the global selection of Torsvik et al. (2012), who also provide the parameters for rotating the poles into Australian coordinates (Figures 7-C; Supplementary Table  S3). We selected this specific 10 Myr time-window because it is approximately centered on the age of the Sommet-Khian section (∼51-46 Ma). The resulting paleolatitude (33.4°S) is significantly lower than the 45.8°S calculated using the dataset from Sommet-Khian. The same applies for the paleolatitude of ∼38°S predicted by the dataset of Veevers and Li (1991). This implies that either the tectonic strain affecting the Sommet-Khian sediments pervasively affected the dataset (even if it seems acceptable within the statistically determined bounds) or the lack of autochthonous (i.e., Australian) data for the Eocene results in an erroneous paleolatitude estimation for the area. We are more inclined to consider that the problem resides in the Sommet-Khian directions, because the AMS data clearly indicate a nonsedimentary fabric. We stress the fact that without AMS analysis, the dataset could have been considered reliable, because of the presence of several reversals and a positive reversal test. A collection of autochthonous Eocene data from either continental Zealandia or Australia might add new constraints to resolve the paleogeography of the southwest Pacific area. Our approach highlights that AMS measurement should be always presented to help assessing the quality of the data.

CONCLUSION
Finite strain in sedimentary rocks induced by tectonism can deviate primary paleomagnetic directions. Using different simulated deformations on distributions drawn according to the TK03. GAD paleosecular variation model, we show how finite strain affects the shape of the directions distribution, which can be used as an indicator of the reliability of paleomagnetic directions. Only if the shape agrees with the expected one based on the TK03.GAD model, a direction distribution can be regarded as of high quality. Otherwise, a "flawed" distribution shape is indicative of post depositional changes requiring further investigation. We apply this concept to a recently published early-middle Eocene paleomagnetic and AMS dataset from New Caledonia. Despite the high quality of the directions and a positive reversal test, the paleomagnetic inclination from New Caledonia is considerably steeper than the one expected from published reference datasets, which are however limited in number and mostly rotated into Australian coordinate from other plates. Paleomagnetic directions from New Caledonia are likely affected by some amount of finite strain, as revealed by a weak prolate AMS fabric of tectonic origin. The distribution shape appears to be in agreement with the TK03.GAD model within a (large) 90%-bootstrap-determined confidence boundaries. Nevertheless, the presence of an AMS tectonic fabric is enough to undermine the reliability of the data. New Eocene data from continental Australia and/or Zealandia should solve the conflict between the autochthonous paleomagnetic dataset and the reference data, which, as is, results in a considerable paleolatitude discrepancy. More importantly, we suggest that AMS measurements should be always measured and presented together with directional datasets when performing paleogeographic reconstructions based on paleomagnetic data from sedimentary rocks.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.