Evolution of part of the <plate> - <plate> boundary in the NNR reference frame

by <your name>

Code by Vince Cronin, revised 4 April 2017

Parts of this code that you will need to change are highlighted in yellow.  You will also need to add appropriate comments and explanations, as indicated by the text enclosed in angle brackets <like this> throughout the notebook, which are not highlighted in yellow.  Delete this cell from your final version of this code that you submit for credit.

1  Introduction

<Explain the purpose of this notebook, what plate boundary section you are working on, which two plates are associated with this boundary, and why you chose to work on the particular plate boundary.  Why is it interesting to you, and why might it be geologically/tectonically interesting?>

<text>

2  Input data

2.1  Angular velocity data for plates relative to NNR reference frame

2.1.1  Instantaneous plate motion in an NNR reference frame from Argus et al. (2011)

<Describe the following data from Argus and others (2011) indicating instantaneous plate motion in an NNR reference frame>

nnr=(({{Antarctica, an, 65.42, -118.11, 0.250}, {Australia, au, 33.86, 37.94, 0.632}, {Caribbean, ca, 35.20, -92.62, 0.286}, {Cocos, co, 26.93, -124.31, 1.198}, {Eurasia, eu, 48.85, -106.50, 0.223}, {NorthAmerica, na, -4.85, -80.64, 0.209}, {Nubia, nb, 47.68, -68.44, 0.292}, {Nazca, nz, 46.23, -101.06, 0.696}, {Pacific, pa, -63.58, 114.70, 0.651}, {SouthAmerica, sa, -22.62, -112.83, 0.109}, {Somalia, sm, 49.95, -84.52, 0.339}}));

2.1.2  Plates and plate boundary segment being studied

<text>

Figure 1.  Boundary segment between <insert plate name> and <insert plate name> from GeoMapApp (http://www.geomapapp.org), showing the approximate current location of the boundary (black polyline).  <insert the rest of the caption here>

In the code below, one plate is designated as plate A and the other is plate B.  (This allows the basic code to be used for different plate pairs without extensive revision.)

Please identify your choice of plate A:  <insert name here -- example: Australia>.

Please identify your choice of plate B:  <insert name here -- example: Somalia>.

2.1.3  Define the full angular velocity vector eΩplate for each plate relative to an external reference frame (NNR)

<Explain the following...>

eΩA=nnr[[2,5]]*convert2Cart[nnr[[2,3]],nnr[[2,4]]];

<Explain the following...>

eΩB=nnr[[11,5]]*convert2Cart[nnr[[11,3]],nnr[[11,4]]];

2.2  Location data along the boundary

<Reference points along the current plate boundary (Bird, 2003;  GeoMapApp).  In any pair of coordinates, the first one is the latitude and the second is the longitude -- north latitude and east longitude are positive; south latitude and west longitude are negative.  Some example data for the ridge between the Australian and Somalian Plates are provided.>

inputBoundary={{-14.741,66.375},{-13.949,65.883},{-13.343,66.718},{-12.769,66.313},{-12.151,65.610},{-11.618,66.357}};

2.3  Model times

This code is written to produce graphical and numeric answers for evenly spaced model times.  So, for example, if you want to model a boundary in 5 Myr steps, you could begin at -5 Myr (1 step in the past) and go to +10 Myr (two steps in the future).  One of the model times must be 0 Myr.

timeStep=5;

pastTimeSteps=1;

futureTimeSteps=1;

2.4  Choices that affect the graphic output

2.4.1  Define the length of the north arrow

<Explain the following...>

northArrowScale=0.1;

2.4.2  Choose the distance from Earth’s center from which the boundary is viewed

Variable viewDistance is the distance to the point from which the boundary is viewed looking toward the center of Earth, measured in earth radii.  The value assigned to this variable must be greater than 1.

viewDistance=1.1;

3  User-defined functions

<Supply a sentence or so of introductory text that answers the following question:  What is this section for? >

The user-defined function convert2Cart takes the latitude (lat) and longitude (long) of a point on Earth’s surface and computes the unit location vector to that point in a reference frame whose origin is Earth’s center, the intersection of the equator and the prime meridian (0° latitude, 0° longitude) is along the positive X axis, the positive Y axis passes through 0° latitude and 90°E longitude, and the north pole is along the positive Z axis.

convert2Cart[lat_, long_] := {Cos[lat Degree] Cos[long Degree], Cos[lat Degree] Sin[long Degree], Sin[lat Degree]};

<Supply text that answers the following question:  What does the following user-defined function do?>

unitVect3D[vect_]:={(vect[[1]]/Norm[vect]),(vect[[2]]/Norm[vect]),(vect[[3]]/Norm[vect])};

<Supply text that answers the following question:  What does the following user-defined function do?>

findGeogCoord[vect_]:=Module[{lat,long,a,b,c,d,e,f},a=ArcSin[vect[[3]]];b={vect[[1]],vect[[2]],0};c=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),{1,1,0},{vect[[1]]/Norm[b],vect[[2]]/Norm[b],0}];d={1,0,0};e=VectorAngle[c,d];f=If[(vect[[2]]<0),(-e),(e)];lat=a(180/π);long=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),0,(f(180/π))];{lat,long}];

The user-defined function northArrow takes the average of all the unit location vectors of points along the plate boundary (viewPoint) and a user-defined scale (nArrowScale), generally ranging from about 0.05 to 0.1, and finds a suitable endpoint for a vector pointing toward north in the output graphics.

northArrow[viewPoint_,nArrowScale_]:=Module[{north,center,normal,northward,answer},north={0,0,1};center={0,0,0};normal=Cross[north,viewPoint];northward=Cross[viewPoint,normal]*nArrowScale;answer={center,northward};answer];

<Supply text that answers the following question:  What does the following user-defined function do?>

circMotion[x_,angVelVect_,dT_]:=Module[{north,rotPole,w,xPole,yPole,m1,m2,m3,answer},north={0,0,1};rotPole=unitVect3D[angVelVect];w=Norm[angVelVect];xPole=unitVect3D[rotPole\[Cross]north]; yPole=unitVect3D[rotPole\[Cross]xPole];m1={{xPole[[1]],xPole[[2]],xPole[[3]]},{yPole[[1]],yPole[[2]],yPole[[3]]},{rotPole[[1]],rotPole[[2]],rotPole[[3]]}};m2={{Cos[(w dT)Degree],-Sin[(w dT)Degree],0},{Sin[(w dT)Degree],Cos[(w dT)Degree],0},{0,0,1}};m3=Inverse[m1];answer=m3.m2.m1.x;answer];

The user-defined function findPtOnNewBoundary takes the unit location vector of a displaced point on one plate (ptA), the unit location vector of the corresponding displaced point on the other plate (ptB), and the unit location vector to the pole of instantaneous relative motion between the two plates (cartRelativePole) and locates the point that is half way between ptA and ptB and whose angular distance from cartRelativePole is the average of the angular distances from cartRelativePole to ptA and ptB.

findPtOnNewBoundary[ptA_,ptB_,cartRelativePole_]:=Module[{geogX,geogY,geogZ,avRefPt,newX,newY,newZ,rotAngle,tM,rotZ,invTM,answer},geogX={1,0,0};geogY={0,1,0};geogZ={0,0,1};avRefPt=unitVect3D[{(ptA[[1]]+ptB[[1]])/2,(ptA[[2]]+ptB[[2]])/2,(ptA[[3]]+ptB[[3]])/2}];newZ=unitVect3D[Cross[cartRelativePole,avRefPt]];newX=cartRelativePole;newY=unitVect3D[Cross[newZ,newX]];rotAngle=(VectorAngle[ptA,newX]+VectorAngle[ptB,newX])/2;tM=(({{Dot[newX,geogX], Dot[newX,geogY], Dot[newX,geogZ]}, {Dot[newY,geogX], Dot[newY,geogY], Dot[newY,geogZ]}, {Dot[newZ,geogX], Dot[newZ,geogY], Dot[newZ,geogZ]}}));invTM=Inverse[tM];rotZ=(({{Cos[rotAngle], -Sin[rotAngle], 0}, {Sin[rotAngle], Cos[rotAngle], 0}, {0, 0, 1}}));answer=invTM.rotZ.tM.newX;answer];

4  Computation

4.1  Convert geographic location data along plate boundary to unit location vectors

Variable refPts is the set of all points along the boundary of interest, expressed as unit location vectors.  Each point is converted from the input geographic coordinates to Cartesian coordinates using the function convert2Cart.

refPts=Table[convert2Cart[inputBoundary[[i,1]],inputBoundary[[i,2]]],{i,1,Length[inputBoundary]}];

4.3  Define the full angular velocity vector observing.plateΩmoving.plate for one plate relative to the other plate

<Explain the following...>

plateAΩplateB=eΩB-eΩA;

<Explain the following...>

plateAωplateB=Norm[plateAΩplateB];

<Explain the following...>

cartRelPole={plateAΩplateB[[1]]/plateAωplateB,plateAΩplateB[[2]]/plateAωplateB,plateAΩplateB[[3]]/plateAωplateB};

<Explain the following...>

geogRelPole=findGeogCoord[cartRelPole];

4.4  Model time

<Explain the following...>

pastTime=0-(pastTimeSteps*timeStep);

futureTime=futureTimeSteps*timeStep;

currentTime=0;

<Explain the following...>

modelTime={pastTime,currentTime,futureTime};

4.5  Compute boundary locations

<Explain the following...>

4.5.1  Compute location vectors for the reconstruction of the past boundary

<Explain the following...>

Variable pastPlateA is a list of points derived by taking the current plate boundary and rotating it around the pole of rotation for that plate by an angle that results from multiplying pastTime by the angular speed.

pastPlateA = Table[circMotion[refPts[[i]], eΩA, pastTime],{i, 1, Length[refPts]}];

Variable pastPlateB is a list of points derived by taking the current plate boundary and rotating it around the pole of rotation for that plate by an angle that results from multiplying pastTime by the angular speed.

pastPlateB = Table[circMotion[refPts[[i]], eΩB, pastTime],{i, 1, Length[refPts]}];

Variable pastBoundary is a list of points that represents the plate boundary at pastTime, and that is (1) half way between pastPlateA and pastPlateB, and (2) the average distance from the pole of relative motion between the two plates and both pastPlateA and pastPlateB.

pastBoundary=Table[findPtOnNewBoundary[pastPlateA[[i]],pastPlateB[[i]],cartRelPole],{i,1,Length[refPts]}];

4.5.2  Compute location vectors for the reconstruction of the present boundary in light of its evolution over the preceding timeStep Myr

<Explain the following...>

outPlateA1 = Table[circMotion[pastBoundary[[i]], eΩA, timeStep],{i, 1, Length[pastBoundary]}];

outPlateB1 = Table[circMotion[pastBoundary[[i]], eΩB, timeStep],{i, 1, Length[pastBoundary]}];

currentBoundary=Table[findPtOnNewBoundary[outPlateA1[[i]],outPlateB1[[i]],cartRelPole],{i,1,Length[refPts]}];

comboPlateA1={currentBoundary,outPlateA1};

comboPlateB1={currentBoundary,outPlateB1};

4.5.1  Compute location vectors for the reconstruction of the present boundary in light of its evolution over the preceding two timeStep Myr

<Explain the following...>

outPlateA2 = Table[circMotion[comboPlateA1[[i,j]], eΩA, timeStep],{i, 1, Length[comboPlateA1]},{j,1,Length[refPts]}];

outPlateB2 = Table[circMotion[comboPlateB1[[i,j]], eΩB, timeStep],{i, 1, Length[comboPlateB1]},{j,1,Length[refPts]}];

futureBoundary=Table[findPtOnNewBoundary[outPlateA2[[1,i]],outPlateB2[[1,i]],cartRelPole],{i,1,Length[refPts]}];

4.7  Prepare the graphics files

4.7.1  Define the view vector

Variable averageVV is the point from which the boundary is viewed looking toward the center of Earth, as plotted using the built-in Mathematica function Graphics3D

averageVV=viewDistance*(Sum[refPts[[i]],{i,Length[refPts]}]/Norm[Sum[refPts[[i]],{i,Length[refPts]}]]);

4.7.2  Compile the graphics files for past, present, and future boundaries

Note:  Some useful (easily visible) built-in colors for use in the graphics include the following:  Red, Green, Blue, Black, Gray, Cyan, Magenta, Brown, Orange, Pink, and Purple.

<Explain the following...>

outPast = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[pastBoundary[[i]],{i,1,Length[refPts]}]],Black,Dashed,Line[Table[refPts[[i]],{i,1,Length[refPts]}]],Red,Line[Table[pastPlateA[[i]],{i,1,Length[refPts]}]],Line[Table[pastPlateB[[i]],{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

<Explain the following...>

outPresent = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[currentBoundary[[i]],{i,1,Length[currentBoundary]}]],Orange,Line[Table[outPlateA1[[i]],{i,1,Length[refPts]}]],Line[Table[outPlateB1[[i]],{i,1,Length[refPts]}]],Black,Dashed,Line[Table[pastBoundary[[i]],{i,1,Length[pastBoundary]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

<Explain the following...>

outFuture = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[futureBoundary[[i]],{i,1,Length[futureBoundary]}]],Orange,Line[Table[outPlateA2[[2,i]],{i,1,Length[refPts]}]],Line[Table[outPlateB2[[2,i]],{i,1,Length[refPts]}]],Blue,Line[Table[outPlateA2[[1,i]],{i,1,Length[refPts]}]],Line[Table[outPlateB2[[1,i]],{i,1,Length[refPts]}]],Black,Dashed,Line[Table[pastBoundary[[i]],{i,1,Length[refPts]}]],Line[Table[currentBoundary[[i]],{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

5  Output

<Explain...>

Show[outPast]

Figure 2.  Boundary modeled at <insert past time> Myr (solid black polyline).  The north arrow is near the current location of the boundary (dashed black polyline).  The dashed red polylines show the current edges of the two plates rotated back by -<startTime> Myr, so the actual plate boundary at that time was between these two dashed polylines. <insert the rest of the caption here>

Show[outPresent]

Figure 3.  Boundary as it is today (solid black polyline).  Plate boundaries from <startTime> Myr ago are shown as orange polylines.  The dashed black polyline is the modeled location of the plate boundary at <startTime> Myr ago.  <insert the rest of the caption here>

Show[outFuture]

Figure 4.  Boundary modeled at <insert future time> Myr (solid black polyline).  The north arrow is near the current location of the boundary (dashed black polyline).  The blue solid polylines show the edge of the <timeStep> Myr crust (as age is measured at <futureTime> Myr), and the orange polylines show the edge of the <two timeSteps> Myr crust.  <insert the rest of the caption here>

Check the relative pole and speed between the two plates.

Compare the computed latitude, longitude, and angular speed of these plates relative to each other (given below) with those published in DeMets et al (2010).

latRelPole=geogRelPole[[1]]

longRelPole=geogRelPole[[2]]

angularSpeed=plateAωplateB

6  Discussion, Conclusions

<text>

7  References and data sources

Argus, D.F., Gordon, R.G., and DeMets, C., 2011, Geologically current motion of 56 plates relative to the no-net-rotation reference frame:  Geochemistry Geophysics Geosystems, v. 12, Q 1101, 13 p., doi:  10.1029/2011GC003751.

Bird, P., 2003, An updated digital model of plate boundaries:  Geochemistry Geophysics Geosystems, v. 4, no. 3, p. 1027, doi:10.1029/2001GC000252.  Available online via http://peterbird.name/publications/2003_PB2002/2003_PB2002.htm

Cronin, V.S., 2017, A primer on the kinematics of the lithosphere and continental crust [draft]:  accessible on the web via http://CroninProjects.org/Vince/PlateKinematics/KinematicsPrimer/index.htm

Cronin, V.S., 1991, The cycloid relative-motion model and the kinematics of transform faulting, in Hilde, T.W.C., and Carlson, R.L., [editors], Proceedings of the 1987 Geodynamics Symposium: Tectonophysics, v. 187, p. 215-249; corrections to publisher errors, 1991, Tectonophysics, v. 192, no. 3/4, p. 401.

Cronin, V.S., 1988, Cycloid tectonics -- A kinematic model of finite relative plate motion:  College Station, Texas A&M University, Ph.D. dissertation, 118 p.

Cronin, V.S., 1987, Cycloid kinematics of relative plate motion: Geology, v. 15, p. 1006-1009.

DeMets, D. C., Gordon, R.G., and Argus, D.F., 2010, Geologically current plate motions: Geophysical Journal International, v. 181, p. 1–80, doi:10.1111/j.1365-246X.2010.04491.x.

Some web resources to investigate:

http://geoscience.wisc.edu/~chuck/MORVEL/

http://www.virtualocean.org/

http://www.geomapapp.org/

http://worldwind.arc.nasa.gov/java/