Cycloid finite motion of points on Carribean plate relative to North American plate

by Joseph Thangraj

Code by Vince Cronin, revised 12 April 2011

1  Introduction

This notebook evaluates and plots the motion of one plate as observed from the other plate, when both plates move about different poles which don’t move relative to each other.The Carribean plate is surrounded by Panama, North American, South American, North Andes and Cocos plates. It has a right lateral transform boundary in the north with the North American Plate, left lateral transform boundary in the south with the South American Plate. The Carribean plate acts as the subducting plate with Panama plate and the Norht Andean Plate, while the Cocos plate subducts under the Carribean Plate. These complex plate motions make the Carribean plate really interesting to study.

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)

Angular velocities of the plates relative to the No-Net-Rotation 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

Figure 1.  Boundary segment between North America and Carribean from GeoMapApp (http://www.geomapapp.org), showing the approximate current location of the boundary (black polyline).

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.)  Plate A is the frame of reference from which the motion of plate B is observed.

Please identify your choice of plate A:  North American plate

Please identify your choice of plate B:  Carribean plate

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

Multiplies the angular velocity with cartesian co-ordinates of the point at plate A.

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

Multiplies the angular velocity with cartesian co-ordinates of the point at plate B.

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

2.2  Location data along the boundary

Location data along the Carribean and North American Plate Boundary.

inputBoundary={{14.877,-90.158},{15.403,-88.893},{16.619,-86.396},{17.292,-84.445},{17.661,-82.477},{17.895,-81.773},{18.48,-81.756},{19.012,-81.738},{19.029,-81.088},{19.178,-79.717},{19.576,-78.047},{20.006,-74.988},{20.237,-72.791}};

2.3  Model times

This code is written to produce graphical and numeric answers for evenly spaced model times, beginning at 0 Myr (today) and continuing in the future by some integral number of time steps.  So, for example, if you want to model a boundary in 5 Myr steps, you could begin at 0 Myr and go to +20 Myr (four steps in the future).

timeStep=5;

futureTimeSteps=3;

2.4  Choices that affect the graphic output

2.4.1  Define the length of the north arrow

Set the size of the north arrrow

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

This section defines all the user-defined functions for use later.

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]};

Normalizes the vector, hence creating a unit vector

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

Gets latitude and longitudes from input cartesian co-ordinates.

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];

Cycloid calculates the motion of the point relative to two poles that are not moving relative to each other. It takes velocity of the observer’s plate, velocity of the moving plate and time as input.

cycloid[(*list of global variables*)(*obsΩ is the total velocity vector of the observer's plate in the external reference frame*)obsΩ_,(*movΩ is the total velocity vector of the moving plate in the external reference frame*)movΩ_,(*refPt is the unit location vector to the reference point whose finite motion is to be modeled*)refPt_,(*t is time in Myr*)t_]:=Module[{(*list of local variables*)a,b,c,d,e,f,g,h0,h1,h,i0,i1,i,j0,j1,j,k,l,m1, m2, m3, m4, m5,n,p,answer},(*geographic X axis -- equator at prime meridian*)a={1,0,0};(*geographic Y axis*)b={0,1,0};(*geographic Z axis -- north pole*)c={0,0,1};(*eωobs*)d=Norm[obsΩ];(*pole eΩobs -- observer coordinate system (CS) Z axis*)e={(obsΩ[[1]]/d),(obsΩ[[2]]/d),(obsΩ[[3]]/d)};(*eωmov*)f=Norm[movΩ];(*pole eΩmov -- other plate CS Z axis*)g={(movΩ[[1]]/f),(movΩ[[2]]/f),(movΩ[[3]]/f)};(*vector coincident with observer CS Y axis*)h0=e\[Cross]g;(*length of h0*)h1=Norm[h0];
(*unit vector along observer CS Y axis*)h={(h0[[1]]/h1),(h0[[2]]/h1),(h0[[3]]/h1)};(*vector coincident with observer CS X axis*)i0=h\[Cross]e;(*length of i0*)i1=Norm[i0];
(*observer CS X axis*)i={(i0[[1]]/i1),(i0[[2]]/i1),(i0[[3]]/i1)};(*vector coincident with other plate CS X axis*)j0=h\[Cross]g;(*length of j0*)j1=Norm[j0];
(*other plate CS X axis*)j={(j0[[1]]/j1),(j0[[2]]/j1),(j0[[3]]/j1)};(*other plate CS Y axis -- same as observer CS Y axis*)k=h;(*m1 is the 3x3 transformation matrix from the geographic CS to the other plate CS*)m1={{j[[1]],j[[2]],j[[3]]},{k[[1]],k[[2]],k[[3]]},{g[[1]],g[[2]],g[[3]]}};(*m2 is a 3x3 Z-rotation matrix for a rotation of (-eωmov*t) degrees around pole eΩmov*)m2={{Cos[(-f Degree) t],Sin[(-f Degree) t],0},{-Sin[(-f Degree) t],Cos[(-f Degree) t],0},{0,0,1}};(*k is the angle between d and j*)l=VectorAngle[e,j];(*n is a variable that helps to differentiate between supplimentary angles*)n=If[l<(π/2),1,-1];(*p is the angle between f and d that has the proper sign to be used in matrix m3*)p=(n)(VectorAngle[g,e]);(*m3 is a 3x3 Y-rotation matrix for a rotation of p radians around the common Y axis of the observer CS and the other plate CS*)m3={{Cos[p],0,-Sin[p]},{0,1,0},{Sin[p],0,Cos[p]}};(*m4 is a 3x3 Z-rotation matrix for a rotation of (eωobs*t) degrees around pole eΩobs*)m4={{Cos[(d Degree) t],Sin[(d Degree) t],0},{-Sin[(d Degree) t],Cos[(d Degree) t],0},{0,0,1}};(*m5 is the 3x3 transformation matrix from the observer CS to the geographic CS*)m5={{i[[1]],h[[1]],e[[1]]},{i[[2]],h[[2]],e[[2]]},{i[[3]],h[[3]],e[[3]]}};(*answer is the product of multiplying the input location vector to a point on the other plate (refPt) by matrices m1, m2, m3, m4, and m5, resulting in the location vector to the point in "t" Myr in a version of the geographic CS that is fixed to the observer's plate*)answer=m5.m4.m3.m2.m1.refPt;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

Angular velocity of plate A relative to plate B

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

Calculate the norm of the angular velocity of plate A relative to plate B

plateAωplateB=Norm[plateAΩplateB];

Calculate unit angular velocity vector

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

Find latitude and longitude from the angular velocity vector

geogRelPole=findGeogCoord[cartRelPole];

4.4  Model time

Set future time from time step and future time steps

futureTime=futureTimeSteps*timeStep;

Create a table of future time steps

modelTime=Table[i,{i,0,futureTime,timeStep}];

4.5  Compute boundary locations

Calculates the cycloid motion of points on the boundary

cycOut = Table[cycloid[eΩA, eΩB, refPts[[i]], modelTime[[j]]], {j,1,Length[modelTime]},{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 file

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.

boundaryOut plots the boundary with time and trajectory plots the trajectory of the points on the boundary with time.

boundaryOut = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[cycOut[[1,i]],{i,1,Length[refPts]}]],Red, Line[Table[cycOut[[j,i]],{j,2,Length[cycOut]},{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

trajectoryOut = Graphics3D[{Arrow[northArrow[averageVV,northArrowScale]],Black,Line[Table[cycOut[[1,i]],{i,1,Length[refPts]}]],Orange, Line[Table[cycOut[[j,i]],{i,1,Length[refPts]},{j,1,Length[cycOut]}]],Blue,Dashed,Line[Table[cycOut[[Length[cycOut],i]],{i,1,Length[refPts]}]]}, AspectRatio->Automatic, Boxed->False,ViewVector->{averageVV,{0,0,0}}];

5  Output

Plot of the Carribean plate boundary as seen from the North American Plate Boundary

Show[boundaryOut]

Figure 2.   Present boundary (solid black polyline) and modeled location of those same points on the Caribbean plate as observed from the North American plate at 5 Myr intervals into the future (solid red polylines).    Black arrow indicates direction toward north.

Plot of the trajectory of the points on the boundary

Show[trajectoryOut]

Figure 3.  Present boundary (solid black polyline) and modeled trajectories of those same points on the Caribbean plate as observed from the North American plate at 5 Myr intervals into the future (solid orange polylines) to their modeled positions at 15 Myr (dashed blue polyline).    Black arrow indicates direction toward north.

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 today (given below) with those published in DeMets et al (2010).

latRelPole=geogRelPole[[1]]

longRelPole=geogRelPole[[2]]

angularSpeed=plateAωplateB

The latitude, longitude and angular speed in DeMets et al (2010) are 73.9,-147.4 and 0.190

6  Discussion, Conclusions

The results of the boundary as observed from the North American plate confirm that the boundary is a left-lateral strike slip boundary. From the plot of the trajectory of the points, we can see that the strike of the transform boundsary has been changing with time. The change in strike is  in the anticlockwise direction. This change in strike will widen the trough between the two plates.

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/