Working with Paleomagnetic Data in GPlates

Authors: Trond H. Torsvik1,2 & Robin J. Watson2

1Physics of Geological Processes & Geosciences, University of Oslo, Norway

2Centre for Geodynamics, NGU, Trondheim, Norway

Working with Paleomagnetic Data in GPlates

Background

Tutorial Files

Exercise 1: Reconstruct a continent according to a palaeomagnetic pole

Exercise 2: Reconstruct multiple continents with palaeomagnetic poles

Exercise 3: Interactively modify a reconstruction: Updating the rotation file

Exercise 4: Reconstructing ‘Laurussia’ from a combined APW path

Exercise 5: The making of Pangea

Exercise 6: The Plume Generation Zone method

References

Background

Figure 1. Inclination of the Earth’s magnetic field and how it varies with latitude

The Earth’s magnetic field is described by its inclination (angle with respect to the local horizontal plane), declination (angle with respect to the Greenwich meridian) and field-strength. The inclination of the Earth's field varies systematically with latitude (Fig. 1), which is of prime importance for palaeomagnetic reconstructions. At the north magnetic pole the inclination of the field is +90° (straight down), at the Equator the field inclination is zero (horizontal) pointing north and at the south magnetic pole the inclination is -90° (straight up). The magnetic north and south poles currently differ from the geographic north and south poles by 11.5° because the magnetic axis is inclined from the geographic (= rotation) axis. The magnetic axis, however, is slowly rotating/precessing around the geographic axis (this is known as secular variation), and over a period of a few thousand years it is hypothesised that the averaged magnetic poles correspond reasonably with the geographic poles. This is known as the geocentric axial dipole (GAD) hypothesis.

Figure 2. Reconstruction of Baltica with paleomagnetic poles. Inset reconstruction uses the south Paleomagnetic pole

Based on measurement of the remanent inclination we can calculate the ancient latitude for a continent when the rock formed from the dipole formula tan(I) = 2 * tan(latitude). In addition, the remanent declination, which deviates from 0° or 180° (depending of the polarity of the Earth’s magnetic field), provides information about the rotation of a continent. The inclination and declination change with the position of the sampled rock on the globe, but the position of the magnetic pole of a geocentric axial dipole is independent of the locality where the rock acquired its magnetisation. Thus, it is practical to calculate pole positions in order to compare results from various sites or to perform plate-tectonic reconstructions.

Ideally, as a time average, a palaeomagnetic pole (calculated from declination, inclination and the geographic site location - see Box 1) for a newly formed rock will correspond with the geographic north or south pole. If a continent moves later, the palaeomagnetic pole must move with the continent. To perform a reconstruction with palaeomagnetic poles we therefore have to calculate the rotation (Euler) pole and angle which will bring the palaeomagnetic pole back to the geographic north or south pole, and then rotate the continent by the same amount. In our example (Fig. 2), a palaeomagnetic pole (latitude=49.3°N, longitude=152.3°E), calculated from the situation depicted in the diagram, will position the Baltica palaeo-continent (most of northern Europe eastward to the Urals) at latitudes between 15-50°N, causing the city of Oslo to have been located at 24°N in the Late Permian. Because the current latitude of Oslo is 60°N, Baltica must have drifted northwards since the Permian.

Palaeomagnetic data can only constrain latitude (based on inclination) and the amount of angular rotation (based on declination). Because the palaeolongitude is unknown, we can position a continent at any longitude we wish subject to other geological constraints. However, this degree of freedom can be minimized by selecting an appropriate reference plate; in other words, if one can determine which plate (or continent) has moved longitudinally the least since the time represented by a reconstruction, then that plate should be used as the reference plate (Burke & Torsvik 2004; Torsvik et al. 2008). Africa is the most appropriate candidate to minimize longitude uncertainty, and thus in the GPlates rotation file, all plates are first reconstructed relative to Africa with appropriate plate-circuits, and then a global apparent polar wander (APW) path constructed in African co-ordinates is used to position all plates at their correct latitude and best possible longitude.

In addition to longitude, we cannot tell in old rocks whether a palaeomagnetic pole is a South or North pole. In Figure 2a we assumed that the pole was a North pole, but if we used a South pole, Baltica would be positioned in the southern hemisphere, and with a geographically inverted orientation (Fig. 2b). Hence, there is freedom to select north or south poles when producing reconstructions, placing the continent in an opposite hemisphere and rotating by 180 degrees — this is not a major problem since Mid-Late Palaeozoic times.

APW paths represent a convenient way of summarising palaeomagnetic data for a continent or terrane instead of producing palaeogeographic maps at each geological period. APW paths represent the apparent motion of the rotation axis relative to the continent depending on whether one plots the movement of the north or south pole. APW paths can therefore be constructed as north or south paths. To construct an APW path, a set of palaeomagnetic poles of varying geologic age are presented in a single diagram, and a synthetic path is fitted to the incrementing poles. There are three common methods for generating APW paths, (1) spherical splines, (2), running mean (sliding-time window) and (3) the small circle method. In our tutorials we include examples of APW paths that are generated with the running mean method averaged over a time-window of 20 Myr.

Tutorial Files

The tutorial files can be downloaded from https://sites.google.com/site/gplatestutorials/, and are organised into separate folders for each exercise.  The files include:

Coastline and other Vector Features:

101_North_America.shp (North American craton)

102_Greenland.shp (Greenland)

302_Balticat.shp (Baltica)

GondwanaSome.shp (South America, South Africa)

Arabia.shp (Arabia)

Iran.shp (Iran)

SouthChina.shp (South China)

ELIP_CHINA_260Ma.shp (Emeishan large ignesous province)

SMEANSLOW1.shp (1% slow contour in SMEAN tomography model)

Palaeomagnetic Mean Poles (GMAP Format; Torsvik & Smethurst 1999; data based on Torsvik et al. 2008a):

Europe2004_RM_Npoles.vgp (Running mean stable Europe 0-330 Ma)

NorthAmerica2004_RM_Npoles.vgp (Running mean North America 0-330 Ma)

Gondwana2010_RM_NPoles.vgp (Running mean Gondwana in Africa frame)

Laurussia2010_RM_NPoles.vgp (Running mean Laurussia in N. America frame)

Rotation Files (Europe and North America based on palaeomagnetic mean poles given above):

Baltica260PM.rot (Euler rotation for Baltica at 260 Ma)

Baltica_Europe.rot (Euler rotations for Baltica/Stable Europe)

North_America.rot (Euler rotations for Baltica/Stable Europe)

Greenland_vs_NorthAmerica.rot (Relative fits between Greenland and North America)

Laurussia_PM.rot (Combined the 3 files above)

Pangea_PM.rot (Euler rotations for some Pangea elements)

Pangea_A_PM.rot (Euler for some elements in Pangea A)

Exercise 1: Reconstruct a continent according to a palaeomagnetic pole

This is a typical task for a palaeomagnetist who after fieldwork and lengthy laboratory work has finally calculated a mean pole from his/her study and now wants to display the resulting reconstruction of the continent where the samples were obtained.

As an example you have studied Permian 260 Ma volcanic rocks in the Oslo Rift (part of Baltica) and obtained a palaeomagnetic pole of ~ 48.1°N and 156.8°E with an error oval A95=4.3°.

1. Go to File → Open Feature Collection (or use the shortcut Ctrl+O) and select file  Tutorial_1\302_Balticat.shp

2. Go to Features → Create VGP (Fig. 3) and type in properties. Type 48.1 and 156.8 for pole latitude and longitude, type 4.3 for A95 etc. (Fig. 4).

3. Click Next and ‘Create a new Feature Collection’ or save VGP pole to an existing Feature Collection (Fig. 5). If you create a new collection you can give it a name later through File-> Manage Feature Collections (shorcut Ctrl+M).

4. Set Time to 260

Figure 3. Step 2 – How to create a virtual geomagnetic pole (VGP)

Figure 4. Step 2 – How to enter the values for your virtual geomagnetic pole

Figure 5. Step 3 - How to turn your VGP into a new feature collection

Figure 6. Baltica and its VGP

On the screen you should now see Baltica (plate ID 302) and the VGP in its present day position (Fig. 6); at this stage we have not yet performed any reconstructions. In order to do that we need to calculate the rotation pole (Euler pole) that will bring the VGP to the present north pole (assuming that it is a north pole).

5. Enable the 'Choose Feature' tool, either by clicking the 'Choose Feature' icon on the Canvas Tools pane (Fig.7, left), or by using the 'Choose Feature' shortcut key (F).

Figure 7. Selecting Baltica's VGP.

6. Select the VGP feature by clicking on the black spot representing the VGP. Ths spot should change colour to white, indicating that it has been selected, and you should see details of the VGP in the 'Current Feature' panel on the right hand side (Fig. 7).

7. Go to Utilities->Calculate Reconstruction Pole (Fig. 8) The fields in this dialog will be pre-filled with values from the VGP which was selected in step 6. Note that it's not necessary to select a VGP before using the 'Calculate Reconstruction Pole' utility; the user can manually enter any desired VGP details in the dialog, whether or not a VGP has been selected.

Figure 8. Calculate Reconstruction Pole dialog.

8. Select “North Pole” and click on Calculate

This should result in a reconstruction pole of latitude 0­° (the latitude should always be zero), longitude 66.8° and angle 41.9°. In future GPlates versions it will be possible to insert this rotation pole into a rotation file by clicking 'Insert Pole in Rotation Model'. Currently we must do this manually using a text editor (e.g. Wordpad). In this tutorial we have already prepared a rotation file (Fig. 9), and this file should be loaded as follows:

7. Go to File → Open Feature Collection and select the file

Tutorial_1\Baltica260PM.rot

Figure 9. Pre-prepared rotation file using reconstructed pole.

After loading the rotation file, Baltica will be reconstructed to its Permian position, and the pole will coincide with the present north pole (Fig. 10).

Figure 10. Baltica reconstructed to its Permian position, with its pole now coinciding with the present day north pole.

Exercise 2: Reconstruct multiple continents with palaeomagnetic poles

1. Go to File → Open Feature Collection and select all files in the Tutorial_2 folder as shown in Figure 10.

Figure 11. Step 1 – How to load multiple features into GPlates.

In this example we load three plates from ESRI shape files (.shp): North America (Plate Id.=101), Greenland (Plate Id.=102) and Baltica (Plate Id.=302). In addition we have opened APW paths for North America (NorthAmerica2004_RM_Npoles.vgp) and Baltica (Europe2004_RM_Npoles.vgp). These are running mean APW paths (20 Myr window) based on raw palaeomagnetic poles listed in Torsvik et al. (2008) with A95 confidence circles. Based on these APW paths we have already calculated Euler poles, and built a rotation file (Laurussia.rot). The rotation file is the ‘engine’ in GPlates and special attention is therefore required.

GPlates is delivered with a default rotation file but you can build your own using a standard text editor (e.g. WordPad). A rotation file follows the standard ‘PLATES’ format and contains 7 columns  (see BOX 2):

Column 1 : Plate ID (3 digit number)

Column 2 : time of reconstruction

Column 3 : Euler latitude (always 0 when calculated from a palaeomagnetic pole)

Column 4 : Euler Longitude

Column 5 : Euler angle

Column 6 : Plate ID of plate to which rotation is relative

(a plate ID of 000 here indicates that the rotations are relative to the spin-axis; this is always used for palaeomagnetic poles)

Column 7 : comment (always started with an exclamation mark)

In our example file we have generated Euler poles for two plates, namely North America (plate ID 101) and Baltica (plate ID 302) based on palaeomagnetic poles. In addition we have included relative fits between Greenland and North America (Gaina 2002, unpublished; listed in Torsvik et al. 2008a); column 6 in the rotation file therefore contains the value 101 to tell the system that the Euler poles are Greenland (102 in column 2) relative to North America (101 in column 6). This rotation file goes back to 330 Myr; you can type any reconstruction age ≤ 330 and GPlates will generate a rotation for that time by interpolate between the rotations in the file.

2. Type 260 in the ’Time’ window and press <RETURN>.

You should now see something like that shown in Fig. 12, with North America and Baltica reconstructed, and the VGPs coinciding with the North pole.

Figure 12. North America and Baltica reconstructed at 260 Ma.

The poles will only coincide with the North pole when they are reconstructed to their VGP age. Press the ‘play’ button to see how this works.  By default, VGPs are visible in a 5 Myr window around their true age. The VGP visibility can be changed through the 'Set VGP visibility' control in the layers menu (Ctrl+L)  (Figs 13 and 14).

Figure 13. The layers window

Figure 14. Setting the VGP visibility

Exercise 3: Interactively modify a reconstruction: Updating the rotation file

1. Go to File → Open Feature Collection and select all files in the Tutorial_3 folder.

2. Type 260 in the ’Time’ window and press <RETURN>

Figure 15. Step 3 – How to choose a feature

The result of this procedure is already described above but we will now modify the relative positions of Baltica and North America/Greenland. These plates were part of Pangea at this time and Baltica was located next to Greenland. However, since Baltica and North America were reconstructed with palaeomagnetic poles, their true longitudes are undetermined, and we must correct this to produce a ‘sensible’ reconstruction.

3. Select the ’Choose Feature’ tool (Fig. 15) (shortcut F) and click on Baltica.

4. Select the ’Modify Reconstruction Pole’ tool (shorcut P).

5. Since palaeomagnetic poles ideally should give you the correct latitude, we only want to adjust the longitude here. We can do this by checking the 'Constrain Latitude' checkbox (in the right hand panel of Fig. 16).  Now, when you click and drag, the selected feature will move along a line of latitude. When the checkbox is not checked, you  can move the plate freely, and you can rotate it by pressing 'Shift' while dragging. If you have two poles with error circles you can statistically move and rotate as long as error circles overlap.

6. Drag the plate until you are happy with the new fit between Greenland and Baltica (Fig. 16) and click ‘Apply’.

7. The ‘Apply Reconstruction Pole Adjustment’ dialog will now appear (Fig. 17). Click ‘OK’ if you are happy with your adjustment.

8. Your new reconstruction (Euler) pole at 260 Ma is now stored only in memory; for permanent storage (Figs. 18 & 19) go to ‘File->Manage Feature Collection’ and save Laurussia_PM.rot with your own filename.

 Figure 16. Step 6 – How to manually manipulate longitudinal positions

Figure 17. Step 7 – How to apply a new reconstruction pole adjustment

Figure 18. Step 8 – How to open the manage feature collections function to save new features

Figure 19. Step 8 – How to permanently save your new rotation files

Exercise 4: Reconstructing ‘Laurussia’ from a combined APW path

1. Go to File → Open Feature Collection and select all files in the Tutorial_4 folder (Fig. 20).

2. Type 260 in the ’Time’ window and press <RETURN>. Laurentia (North America, Greenland) and Baltica will be reconstructed (Fig. 21).

 Figure 20. Step 1 – How to load feature collections

Figure 21. Step 2 – Reconstruction of North America, Baltica and Greenland

This tutorial will explain how you can build a rotation file consisting of palaeomagnetic data and relative plate circuits to position continents in latitude and orientation. The rotation file (Fig. 22)  was built from a global apparent polar wander (APW) path in North American co-ordinates. (The VGPs can be found in the file Laurussia2010_RM_NPoles.vgp).  We then calculated rotation (Euler) poles for each mean pole in 10 Myr intervals. The first lines in the rotation file describe the movement of North America (plate ID 101) relative to 000 (which represents the spin axis), and the rotations can be used back to 330 Ma (around when Pangea formed).  After the North America rotations, relative plate motions between Greenland and North America (accounting for Labrador Sea/Baffin bay opening in the Late Cretaceous and Tertiary) are given.  Finally, the relative positions of Baltica (now part of Europe, plate ID 302) are given with respect to North America (plate ID 101).

Familiarize yourself with the rotations by typing different reconstruction times between 320 and the present day (Fig. 22).  Use the animation controls (Fig. 23) and see how the reconstruction varies with time. The VGPs should coincide with the north pole when the reconstruction time is exactly the same as the VGP age.

By default, mean VGPs are visible in a 5 Myr window around their true age. The VGP visibility can be changed as described in Exercise 2 (Figs. 13 and 14).

Figure 22. Rotation file

Figure 23. How to manipulate animation reconstruction time

Exercise 5: The making of Pangea

1. Go to File → Open Feature Collection and select all files in the Tutorial_5 folder (Fig. 24).

2. Type 260 in the ’Time’ window and press <RETURN>. Parts of Laurussia (North America, Greenland, Baltica) and parts of Gondwana (Africa, South America) will be reconstructed (Fig. 25). Laurussia and Gondwana constituted the major parts of Pangea at this time but obviously our reconstruction does not look like the Pangea supercontinent when reconstructed from Laurussian (Laurussia2010_RM_NPoles.vgp) and Gondwanan

(Gondwana2010_RM_NPoles.vgp) poles — so what is wrong?

 Figure 24. Step 1 – How to open tutorial 5’s feature collection

Figure 25. Step 2 – Reconstruction at 260 Ma. What is wrong with this picture?

Figure 26. How to choose features

The rotation file used in this exercise contains Euler poles calculated from mean APW paths for Laurussia (in North American co-ordinates - plate ID 101) and Gondwana (in South Africa co-ordinates - plate ID 701) along with relative plate circuits to connect the various plates in this example. Since longitude is not constrained by palaeomagnetic data, Gondwana (South Pangea) and Laurussia (North Pangea) will not have their correct relative longitude, so we must move one of them sideways (E-W) to make a proper Pangea reconstruction.

We have already described how to do this in Exercise 3 but here is the procedure again:

· Select the 'Choose Feature' tool (shortcut F) (Fig. 26) and click on South Africa (this is the plate that the palaeomagnetic data in the rotation file relate to).

· Select the 'Modify Reconstruction Pole' tool (shortcut P).

· Click on the 'Constrain Latitude' checkbox and the 'Highlight children' checkbox in the panel to the right.

· Drag South Africa until you are happy with the new fit between Laurussia and Gondwana (Fig. 27) and click 'Apply'.

· The 'Apply Reconstruction Pole Adjustment' dialog will now appear; click 'OK' if you are happy with your adjustment.

· Your new reconstruction Euler pole at 260 Ma is now stored only in memory; for permanent storage go to

'File->Manage Feature Collections' and save Pangea_PM.rot with your own filename.

Figure 27. How to longitudinally adjust Gondwana.

Exercise 6: The Plume Generation Zone method

Burke & Torsvik (2004) and Torsvik et al. (2006, 2008b) showed that the centers of most large igneous provinces (LIPs) of the past 300 Myr, when restored to their eruption sites, lie radially above one or other of two narrow belts centered on the 1% slow shear wave velocity contour of the SMEAN model of Becker & Boschi (2002) at the core-mantle-boundary (CMB). We can use the remarkable correlation between LIPs and CMB heterogeneities to estimate the longitudes for LIPs that erupted on plates that were not part of Pangea. As an example, the China blocks were unconnected to Pangea in late Permian times and have thus no longitudinal constraint. However, the 258 Ma Emeishan LIP in South China has excellent palaeomagnetic data that position it at 4°N. If this LIP erupted above the ~1% slow contour, there are five possible longitudinal locations where the 4°N line of latitude intersects or touches the 1% slow contour. Pangea covered the African large low shear velocity province (LLSVP) options, leaving only the two options related to the Pacific LLSVP.

Figure 28. In late Permian times Pangea was centered above the African LLSVP (SMEAN model). Pangea did not include N and S China, which were located within the Palaeotethys Ocean. Because S China was not part of Pangea its longitudinal relation to South Africa is unknown and palaeomagnetic data allow it to be placed anywhere in palaeolongitude. However, if the 258 Ma Emeishan LIP was erupted from a plume derived from the 1% slow shear wave velocity contour in the lowermost mantle of one or other of the Earth’s two major large low shear velocity provinces (LLSVP's), the most probable position would be along the western edge of the Pacific LLSVP (see Torsvik et al. 2008b)

A reconstruction with Emeishan erupted above the western margin of the Pacific LLSVP at 140ºE (Fig. 28), is in fact the only likely position because the alternative - that Emeishan erupted above the eastern margin of the Pacific LLSVP - would yield unrealistically high plate convergence rates between Eurasia and South China (25-30 cm/yr) between the Late Permian and the Jurassic.

We will now show how this can all be done in GPlates.

1. Go to File → Open Feature Collection and select all files in the Tutorial_6 folder (Fig. 29).

2. Type 258 in the ’Time’ window (the age of the Emeishan LIP) and press <RETURN>. Parts of Laurussia (North America, Greenland, Baltica) and parts of Gondwana (Africa, South America, Arabia, Iran) will be reconstructed (Fig. 30).

Figure 29. Step 1 - How to open Exercise 6’s feature collections

Figure 30. Step 2 – Reconstruction of Pangea with the Emeishan LIP and SMEAN contour files loaded.

The rotation file (which in this case only contains data covering the interval 300 to 250 Myr) is constructed in such a way that we now see a Pangea A reconstruction (with many pieces not shown) with South China sitting at some distance away from Pangea in the Palaeotethys. All plates are relative to South Africa (plate ID 701) except South China (plate ID 602) which is based on its own palaeomagnetic data. The yellow line is the 1% slow contour in the SMEAN model (this line is fixed relative to the spin-axis and thus will not move at any time).

3. Drag the Globe ('Drag Globe' tool, shortcut D) so that South China appears roughly in the centre and so that you can see the margin of the Pacific LLSVP in yellow. The Emeishan LIP related volcanics are shown as numerous small orange polygons (Fig. 31).

4. Select the 'Choose Feature' tool (Fig. 26) and click on South China.

5. Select 'Modify Reconstruction Pole'.

6. Click on the 'Constrain Latitude' checkbox in the panel to the right

7. Drag South China until you are happy that the Emeishan LIP is centered around the 1% slow SMEAN contour (yellow line in Fig. 32) and click 'Apply'.

8. The 'Apply Reconstruction Pole Adjustment' dialog will now appear; click 'OK' if you are happy with your choice.

9. Your new reconstruction Euler pole at 258 Ma is now stored only in memory. For permanent storage go to

'File->Manage Feature Collections' and save Pangea_A_PM.rot with your own filename.

Figure 31. Step 3 – Readjusting the globe to center it on South China

Figure 32. Step 7 – How to longitudinally adjust South China so it lines up with the SMEAN contour

GPlates also handles raster graphics and we can therefore show the SMEAN model at 2800 km depth together with our shapefile features:

1. Go to File → Import Raster and select tomo-2800.jpg from the Tutorial_6 folder (Fig. 33). Click Open → Next → Next → Finish. See Tutorial 10 for more detail about loading and manipulating rasters in GPlates.

Note, that you may need to change the colouring of loaded features by going to Features → Manage Colouring. Alternatively, you may change the order of the layers (e.g. rasters and loaded coastlines) in the 'Layers' Window by dragging the raster to the bottom of the list of loaded files.

The SMEAN model should now appear (Fig. 34); the two LLSVPs appear as red (hot but dense) regions, and blue regions denote the fastest velocities (slabs?) in the deep mantle.

Figure 33. Step 10 – How to load a tomography raster.

Figure 34. The SMEAN tomography model raster. LLSVP’s correlate well with fast zones.

You may want to make use of one of the map projections to see both of the LLSVPs in the same view (Fig. 35).

Figure 35. Using a Mollweide projection.

References

Becker, T.W., Boschi, L., 2002. A comparison of tomographic and geodynamic mantle models. Geochem. Geophys. Geosyst. 3, 1003, doi:10.1029/2001GC000168.

Burke, K. & Torsvik, T.H. 2004. Derivation of large igneous provinces of the past 200 million years from long-term heterogeneities in the deep mantle. Earth Planet Sci. Lett., 227, 531-538.

Burke, K., Steinberger, B., Torsvik, T.H. & Smethurst, M.A. 2008. Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the Core-Mantle Boundary. Earth and Planetary Sciences, 265, 49-60.

Butler, R.F. 1992. PALEOMAGNETISM: Magnetic Domains to Geologic Terranes. Blackwell Scientific Publications (http://www.geo.arizona.edu/Paleomag/book/).

Torsvik, T.H. & Smethurst, M.A. 1999. Plate Tectonic modeling: Virtual Reality with GMAP. Computer & Geosciences, 25, 395-402.

Torsvik, T.H., Smethurst, M.A., Burke, K. Steinberger, B. 2006. Large Igneous Provinces generated from the margins of the Large Low Velocity Provinces in the deep mantle. Geophysical Journal International, 167, 1447-1460.

Torsvik, T.H., Müller, R.D., Van der Voo, R., Steinberger, B. & Gaina, C. 2008a. Global Plate Motion Frames: Toward a unified model. Reviews of Geophysics, 46, RG3004, doi:10.1029/2007RG000227.

Torsvik, T.H., Steinberger, B., Cocks, L.R.M. & Burke, K. 2008b. Longitude: Linking Earth's ancient surface to its deep interior. Earth Planet Science Letters, 276, 273-283.