Analysis of Simpson Island paleomagnetic data¶

Corresponding Author: Nicholas L. Swanson-Hysell (swanson-hysell@berkeley.edu)

Introduction¶

This IPython notebook contains data analysis on paleomagnetic data developed from ca. 1.1 billion year old lava flows of the Osler Volcanic Group that are part of the North American Midcontinent Rift as well as supporting analyses on other data sets. These analyses accompany a Geochemistry, Geophysics, Geosystems manuscript entitled "Confirmation of progressive plate motion during the Midcontinent Rift's early magmatic stage from the Osler Volcanic Group, Ontario, Canada" by N. L. Swanson-Hysell, A. A. Vaughan, M. R. Mustain and K. Asp. This notebook is part of the supporting online materials and is available at https://github.com/Swanson-Hysell/2014_Swanson-Hysell-et-al_Osler along with the necessary files to execute all of the code. Within this github repository are the following folders:

1. 2014_Osler_Code This folder contains this notebook file as well as necessary libraries.
2. 2014_Osler_Data This folder contains the data generated in the study at the specimen level as well as the flow means that are imported into this notebook for data analysis.
3. 2014_Osler_Manuscript This folder contains the manuscript text and figures.

If you are viewing this supporting material as a PDF document and want to run the code in the notebook you will need to download the code from the github repository and you will also need a Python distribution that includes IPython. There are good instructions for installing IPython here: http://ipython.org/install.html. Alternatively you can view the notebook online with the IPython nbviewer at this link: http://bit.ly/19XRdjJ

Import libraries¶

This code blocks imports necessary libraries that define functions that will be used in the data analysis below. The code below uses the pmag.py and pmagplotlib.py files of the the PmagPy software package (version pmagpy-2.206) authored by Lisa Tauxe (https://github.com/ltauxe/PmagPy). The check_updates and get_version functions of the PmagPy libraries cause errors in the interactive environment of IPython so the code below calls a slightly modified version of those libraries that is included in the Github repository for this work. There are other functions that are necessary for the interactive plotting of paleomagnetic data (some of them edited from PmagPy) that are within the IPmag.py library that is imported in the code block below.

This notebook runs with pylab inline which imports the numpy, scipy and matplotlib functions and allows for the plots to be viewed inline in the IPython notebook (instead of opening up in another window).

Import Simpson Island paleomagnetic data¶

The data file FlowDataAll.csv has data in this format:

SITE, STRAT_HEIGHT, DEC_GEO, INC_GEO, Dec_TC, INC_TC, A95, n, Plat, ABS_Plat, VGP_lat, VGP_long

where:

SITE is the name of the site and corresponds to an individual lava flow where in the first part of the site name (e.g. SI1) corresponds to the name of a stratigraphic section while the second part (e.g. (11.8 to 26.4)) corresponds to the stratigraphic ranges within the measured section over which the flows is exposed.

STRAT_HEIGHT is the cumulative stratigraphic height of the Simpson Island stratigraphy starting with the base of the SI1 section (see map in the main manuscript).

DEC_GEO is the site mean declination in geographic (i.e. in situ) coordinates prior to tilt-correction.

INC_GEO is the site mean inclination in geographic (i.e. in situ) coordinates prior to tilt-correction.

Dec_TC is the site mean declination in tilt-corrected coordinates following correction for bedding tilt.

INC_TC is the site mean inclination in tilt-corrected coordinates following correction for bedding tilt.

A95 is the 95$\%$ Fisher confidence ellipse around the calculated site mean.

n is the number of individually oriented samples that are being included to calculate the site mean.

Plat is the paleolatitude calculated from the site mean using the dipole equation.

ABS_Plat is the absolute value of this paleolatitude.

VGP_lat is the latitude of the virtual geomagnetic pole calculated using the site location and the site mean direction.

VGP_long is the longitude of the virtual geomagnetic pole calculated using the site location and the site mean direction.

These data are displayed in a table below and imported by variable in the cell below the table.

Table of Simpson Island paleomagnetic data¶

The .csv of these data can be downloaded here: SimpsonIsland_OslerData.csv

Comparison between data from the lower 1/3, middle 1/3 and upper 1/3 of Simpson Island stratigraphy¶

Let's consider the data by stratigraphically grouping flow data from the lower third (0 to 1041 meters), middle third (1041 to 2083 meters) and the upper third (2083 to 3124 meters) of the stratigraphy. The Fisher means for these groups are calculated with the resulting dictionaries from the fisher_mean pmag.py fuction being printed below and with the data being displayed on equal area plots with the calculated Fisher means also being plotted.

Test for a common mean between paleomagnetic data from the lower, middle and upper thirds of the stratigraphy¶

The code below will compare the subsets of data using these two statistical tests:

1. The $V_w$ test of Watson (1983) slightly modified from the implementation of PmagPy (Tauxe, 2013). The test is comprised of calculating the $V_w$ statistic for the two data sets and then determining the critical value of $V_w$ through a Monte Carlo simulation. For the simulation, two Fisher-distributed data sets with a common mean are simulated using their precisions ($k_1$ and $k_2$) and the number of points ($N_1$ and $N_2$) of the datasets being evaluated. Then the $V_w$ statistic that comes from comparison between the two simulated data sets is calculated. A large number of simulations are done (default is 1000; it can be set as the NumSims parameter of the function) in order to determine the $V_w$ values that would result by chance through sampling distributions with the same direction. The critical value of $V_w$ is at the 95% level of confidence (i.e. if 1000 $V_w$ values are simulated, it is the 950th one). This implementation of the test also provides the critical angle between the two sample mean directions and the corresponding McFadden and McElhinny (1990) classification.
2. The bootstrap fold test of Tauxe (developed by Tauxe et al., 1991 and implemented per Tauxe (2010)). This approach determines the cumulative distributions of the Cartesian coordinates of bootstrapped means and compares the cumulative distributions to see if the confidence intervals overlap. If they do all overlap, then the two means cannot be distinguished at the 95% level of confidence and they pass the bootstrap test for a common mean. Otherwise, if the two sets of directions are distinct in the X, Y or Z component the two means can be distinguished at the 95% confidence level.

Common mean tests between SI_LowerThird_Directions and SI_MiddleThird_Directions:¶

Common mean tests between SI_MiddleThird_Directions and SI_UpperThird_Directions:¶

Common mean tests between SI_LowerThird_Directions and SI_UpperThird_Directions:¶

Common mean test on directions results summary¶

The above results from common mean statistical tests show that the directions from the lower third cannot be distinguished as being obtained from a distribution distinct from that of the middle third nor can the directions from the upper third be distinguished from those of the middle third. In contrast, directions from the lower third of the stratigraphy can be distinguished from those of upper third at the 95% confidence level.

Common mean tests between SI_LowerThird_Poles and SI_UpperThird_Poles:¶

For completeness, let's also conduct tests for a common mean between the virtual geomagnetic poles (VGPs) from the lower third and upper third of the stratigraphy in addition to the tests on the site mean directions in directional space (declination, inclination).

As when the site declinations and inclinations are considered, the virtual geomagnetic poles from the lower third and upper third of the stratigraphy can be distibuguished at the 95% confidence level according to these tests.

Plotting the data stratigraphically¶

Using the flow means calculated for each portion of the stratigraphy, the paleolatitude and 95% confidence bounds on paleolatitude can be calculated. The stratigraphic bounds (y-axis error bars) and paleolatitude bounds (x-axis error bars) are displayed for points plotted in the middle of the stratigraphic bin and at the mean paleolatitude.

Comparing larger and smaller stratigraphic subsets of the Simpson Island data¶

Comparing data from the lowermost 500 meters to data from the uppermost 500 meters¶

The analysis above demonstrates that the directions from the lower third of the stratigraphy (0 to 1041 meters; 30 flows) and the upper third of the stratigraphy (2083 to 3124 meters; 34 flows) are statistically distinct and thus indicate statistically significant motion to lower latitudes. While it is preferable to have a high total number of flows for such comparisons in order to average out secular variation, let's conduct this same analysis, but for a more restricted range from the bottom 500 meters (17 flows) and top 500 meters (17 flows) of the stratigraphy.

The results of this comparison between the lower 500 meters and upper 500 meters show that they are distinct populations at the 95% confidence level. The difference between the populations is more dramatic in both the Watson V and bootstrap test from the tests on the lower third and upper third subsets that were conducted above. This increase in difference with tighter stratigraphic groups from lower and higher in the sequence is the expected result if the sequence is a record of progressive paleogeographic change.

Comparison between data from the lower 1/2 and upper 1/2 of Simpson Island stratigraphy¶

Calculating and plotting mean pole positions¶

The code below calculates the Fisher mean parameters for poles of the stratigraphically grouped virtual geomagnetic poles (VGPs) and then plots them on a "view from space" globe. Note that while where the textual output says "dec" it actually refers to pole longitude and where it says "inc" it actually refers to pole latitude since the pmag.fisher_mean function is being conducted on VGPs.

Note that in the print out for the mean pole calculations above 'dec' is really the longitude of the pole (Plong) and 'inc' is actually the latitude of the pole (Plat).

Comparison between Halls (1974) data and Simpson Island VGPs¶

Halls (1974) developed paleomagnetic data from the Osler Volcanic Group in the Nipigon Strait region (bibtex citation below; also see http://earthref.org/MAGIC/9518). For the purposes of the previously contributed MagIC database entries, site locations were determined using georeferenced versions of the maps provided in Halls (1974). These data allow for the site means of Halls (1974) to be used to calculate VGPs for purposes of directional statistical tests and for combining into a mean paleopole. The Halls (1974) data are predominantly from flows with reversed magnetization below an angular unconformity on Puff Island that separates the reversedly magnetized flows from younger flows of normal polarity (only a few of which are preserved before the sequence is submerged beneath Lake Superior). The analysis below compares paleomagnetic data from this study with the reversed flows from the Halls (1974) data. The result from this analysis is that the Halls (1974) data is statistically distinct from the lower third of the Simpson Island stratigraphy, but is statistically indistinguishable from the upper third of the Simpson Island stratigraphy. This result fits with stratigraphic considerations that place the sites studied by Halls towards the top of the stratigraphy studied herein at Simpson Island.

The code below imports the directional data from Halls (1974) and than uses this data to calculate virtual geomagnetic poles for each site location.

Calculate VGPs using the Halls (1974) data:

Split the dataframe into one for the normally magnetized sites and one for the reversed polarity sites

Develop a list of unit vectors of the reversed Halls (1974) VGPs in order to conduct common mean tests with the Simpson Island data using pmag.py functions.

Common mean tests between Halls (1974) data and the data from the lower third of the Simpson Island stratigraphy¶

A comparison between VGPs calculated from Halls, 1974 data and data from the lower third of the Simpson Island stratigraphy fails Watson's V and bootstrap tests for a common mean indicating that the populations of directions do not share a common mean at the 95% confidence level.

Common mean tests between Halls (1974) data and the data from the upper third of the Simpson Island stratigraphy¶

A comparison between VGPs calculated from Halls, 1974 data (see http://earthref.org/MAGIC/9518) and data from the upper third of the Simpson Island stratigraphy passes Watson's V and bootstrap tests for a common mean.

Combining the data from the upper third of the Simpson Island stratigraphy and the Halls (1974) reversed polarity data into a single pole for the upper portion of the reversed Osler Group stratigraphy¶

Paleogeographic work outside of this IPython notebook using the software package GPlates.¶

For details on working with paleomagnetic data in GPlates check out this tutorial.

The steps necessary to develop the reconstruction using these Osler Volcanic Group poles are detailed below.

(1) Develop .rot file that rotates Laurentia such that the Osler poles are at geographic north. In the text for the .rot file below the Osler_lowerthird pole is assigned to 1110 Ma, the Osler_middlethird pole is assigned to 1107.5 and the Osler_upperthird pole is assigned to 1105 Ma.

The .rot file has the format:

In this .rot file, the fixed reference frame is 000 and 001 (which could be split into a relative and absolute reference frame that seperates out TPW), while the plate ID for Laurentia is 199.

(2) Open the following feature collections in GPlates: 199-Laurentia nogrid.dat, OslerPoles.gpml and Osler.rot where 199-Laurentia nogrid.dat is the outline of Laurentia in the late Proterozoic and OslerPoles.gpml contains the three calculated Osler paleomagnetic poles. An animation from 1110 to 1105 Ma yields the following reconstruction shown here with a time slice for Osler_lowerthird (red), Osler_middlethird (yellow) and Osler_upperthird (blue).

Plate rate estimate¶

This portion of the analysis calculates plate motion rates and estimates the uncertainty associated with these rates using Monte Carlo simulations that utilize both the uncertainty on the position and age of two paleomagnetic poles. For this analysis we are using these two poles:

1. Pole 1: Osler Volcanic Group (early magmatic stage). The pole is the combined mean of the upper portion of the reversed stratigraphy that uses data from this study and Halls (1974). The date of 1105 $\pm$ 2 Ma is used to constrain the age of the pole since it comes from an extrusive unit in close stratigraphic proximity to the units from which the paleomagnetic data were obtained.
2. Pole 2: North Shore Volcanic Group (main magmatic stage). The North Shore Volcanic Group spans more than 10 million years and subsets of the data need to be considered for a pole to provide precise constraints. The pole used here is calculated from data developed by Tauxe and Kodama (2009) and is the mean of 47 sites from between the 40th Ave icelandite and the Palisade rhyolite within the southwest limb of the North Shore Volcanic Group. The 40th Ave icelandite has a U-Pb date of 1098.4 $\pm$ 1.9 Ma and the Palisade rhyolite has a U-Pb date of 1096.6 $\pm$ 1.7 Ma.

Calculating the NSVG pole¶

The North Shore Volcanic Group (NSVG) is comprised of two main limbs with distinct stratigraphy and radiometric age control. The southwest limb of the NSVG was particularly well-sampled by Tauxe and Kodama (2009) and those sites can be bracketed with age control from the 40th Ave icelandite (Davis and Green, 1997; 1098.4 ± 1.9 Ma) and the Palisade rhyolite (Davis and Green, 1997; 1096.6 ± 1.7 Ma).

Input pole data and simulation sample size¶

The code box below is where the input parameters for the Monte Carlo rate simulation are entered. These parameters are:

• samplesize determines the number of random pole pairs that will be made through Monte Carlo simulations of each pole and its age to determine rate.
• Paleomagnetic pole longitude (plong) and latitude (plat) for each pole
• The A$_{95}$ confidence ellipse (A95), the Fisher precision parameter (kappa), and number of virtual geomagnetic poles used to calculate the pole mean (N).
• The age assigned to the pole and the 1$\sigma$ uncertainty on that age (age_error).

Sampling the age distributions¶

Weighted means calculated from geochronological data are typically assumed to have uncertainity that follows a Gaussian distribution. Using this assumption, we can use the numpy function random to randomly sample a normal distribution. This approach is straightforward for a paleomagnetic pole where we are assigning a single age (e.g. pole1 in this case). It is less straightforward when using two dates that stratigraphically bracket flows from which the paleomagnetic data was obtained without an a priori assumption of the distribution of time between those two dated units. In this analysis, we approximate the age error distribution for pole2 by taking the mean of the dated units that bracket the data as the mean age and then specifying an age uncertainty wherein 2$\sigma$ reaches out to approximately reach the maximum and minimum dates implied by the 2$\sigma$ uncertainties on the older and younger dates respectively. This approach specifies a distribution that is centered on the mean of the two dates while also including the possibility that the actual age of the pole could come from a distribution that spans out to include the older and younger possibilities of each dated unit.

Sampling pole distributions¶

To determine the uncertainty on the pole position and associated paleolatitude, we can take the approach of random sampling from the Fisher distribution with the precision of the VGPs used to develop the pole. To do this in a Monte Carlo fashion, we can grab random samples utilizing the fshdev function from pmag.py which takes random samples from a distribution with a given Fisher precision parameter (kappa) centered on (Dec=0, Inc=90). We can then define a tilt_direction to be the longitude of the pole mean and a tilt_amount to be the conjugate of the latitude of the pole mean and then use the dotilt function of pmag.py to rotate the polar-centered distribution to be centered about the pole mean. Doing this for the same number of VGPs used to calculate the actual pole allows us to then calculate a mean pole. This can then be done repeated (as dictated by the sample size number) to get a population of mean poles that can result from sampling the distribution.

Calculate and graphically display the confidence interval on the minimum rate estimate implied by the Monte Carlo sampled pairs of sampled ages and paleolatitudes from the two poles¶

The pairs of sampled ages and poles can be used to calculate rates (with the total number of rates being set by samplesize in the input box above). The rate calculated above of 24.0 cm/year can now be given confidence bounds by taking 2.5 percentile and 97.5 percentile of the Monte Carlo simulated rates.

We can see here that the 50th percentile from these analysis is the same as the mean rate previously calculated (24.0 cm/yr). The 2.5th and 97.5th percentile give a 95% confidence range of 15.2 to 44.4 cm/yr. The code below generates a plot where 5,000 of the 100,000 (that's what was set as samplesize when code was executed) pole pairs and rates are shown (in A) and a histogram of all the rates is plotted along with the above percentiles being marked. This plot is Figure 4 of the manuscript.

Latitudinal change from early in the early magmatic stage to the main magmatic stage¶

We can use the pole from the lower third of the Simpson Island stratigraphy and the pole calculated above from the NSVG to estimate the magnitude of paleolatitudinal change that occured during the early magmatic stage, the latent magmatic stage and into the main magmatic stage. This paleolatitudinal difference is discussed in the manuscript text.

Lava flow thickness data¶

In the section of the text entitled "Osler Volcanic Group lithologies" the median flow thickness is reported from the subset of flows within measured sections where there is sufficient exposure of the flow base, interior and top to determine thickness. In the text it is reported that the median flow thickness is 4.9 meters with a first quartile thickness of 2.0 meters and third quartile thickness of 9.8 meters. The flow thickness data used in this analysis is presented and analyzed below. Note that not every flow within the measured stratigraphic sections was sampled for paleomagnetic analysis and not every flow on Simpson Island was measured in a stratigraphic section.

The IPmag.py library of functions¶

Some of the functions used in the analysis above come from the imported IPmag.py library that was developed for this data analysis project relying heavily on the pmag.py functions developed by Lisa Tauxe. The code that is within IPmag.py is shown below so that the PDF output of this IPython notebook can document more of the tools used in the analysis without the reader needing to look into the .py files themselves. This library utilizes functions from the pmag.py and pmagplotlib.py libraries of PmagPy version 2.206.