Circular Relative Motion Model

Vince Cronin, 17 April 2010;  revised 20 April 2010

Written in Mathematica, version 7

Introduction

The most important thing to know about the circular relative motion model is that, in general, it does not work for the finite relative motion of plates.  The principal reasons for considering this model are that (a) it is the simplest model of relative motion on a sphere, (b) it has historical significance in the development of plate kinematic ideas, and (c) many people seem to believe that it is useful for finite relative motion.  Consequently, to the minor extent that geoscience textbooks consider modeling of the finite motion of plates at all, textbook authors tend to offer only the circular motion model for students to learn.
There are two reasons to assert that the circular motion model does not work for finite plate motion relative to other plates.  
First, circular finite relative motion requires that the axis of relative motion between the two plates must remain fixed to both plates throughout the time interval.  As soon as we introduce a third plate into the system, we have not one but three axes of relative motion to account for -- each plate associated with two axes.  Unless those two axes are actually coaxial/coincident, they cannot remain fixed relative to the plates whose instantaneous relative motion they describe.  This is called the 3-plate problem of plate tectonics (Cox, 1973).  Given the observed plate system in which none of the axes of instantaneous relative motion appear to be coincident (e.g., Demets and others, 2010), the finite relative motion of plates is not circular.  
Second, the circular motion model would generate oceanic fracture zones that are circular arcs concentric with each other around the axis of relative motion.  In fact, we observe that oceanic fracture zones are not circular, are not concentric, and are not symmetric across the ridge axis.  The shape of oceanic fracture zones appears to be more sigmoid, with the trace on one side of the ridge axis concave in one direction, and the other side concave in the opposite direction.

Assumptions

We assume that Earth is a sphere of unit radius in all but one part of this analysis.  Of course, Earth is not really a sphere.  Earth has an irregular surface with deep ocean basins and tall mountains.  The average radius of Earth is around 6,371.01±0.02  km (Yoder, 1995, p. 8).  The deepest part of the ocean basins, called the Challenger Deep in the Mariana Trench, is ~11.0 km below sea level (that is, ~0.2% of Earth's mean radius), and the summit of Mt. Everest is ~8.848 km above sea level (~0.1% of Earth's mean radius).  Earth's polar radius (6356.7530 km) differs from its mean equatorial radius (6378.1363 km) by only ~21 km, so Earth's shape is different from a sphere by about 1 part in 300 or about 0.3 % (e.g., Cazenave, 1995, p. 36).  
Assuming that Earth is a sphere of unit radius makes the mathematics we use in this problem much simpler than if we used either an ellipsoid or the geoid.  For just one example, the radius of a sphere is the same measured at any point on the sphere's surface, and this is not true for either the reference ellipsoid or the geoid.  Although Earth is a little wider at its equator than along its spin axis, we do not generate serious errors in our calculations by assuming Earth's radius is the same in all directions -- that it is a sphere.   
When we say that Earth is assumed to have unit radius, that means that we set Earth's radius to 1 unit of length.  The one part of this analysis in which we must use the actual size of the Earth in kilometers is in the determination of the magnitude of the tangential velocity vector (i.e., the tangential speed) at a point.

Instantaneous Relative Velocity Models

Geoscientists have generated several models of the instantaneous relative velocities of plates.  By "instantaneous" we mean "over a short time interval."  The available data allow us to define a velocity vector for a given point (that is, the direction and speed of the point), but not the actual path of motion (or finite trajectory) of that point over an extended time interval.
The current standard model for the instantaneous relative motion of plates is the MORVEL model of DeMets and others (2010).  MORVEL is the latest in a series of models developed since the early 1970s, including the models of Chase (1972, 1978), Minster and others (1974), Minster and Jordan (1978), and the NUVEL and NUVEL-1A models of Charles DeMets and others (1990, 1994).  In the MORVEL model (DeMets and others, 2010, table 2), the North American plate moves with an anti-clockwise rotation around a pole located at latitude 48.9°N and longitude 71.7°W with an angular velocity of 0.750°/Myr relative to the Pacific plate.

Input Data

Input only decimal degrees (i.e., not degrees, minutes, seconds).  The sign convention follows:  north latitude is positive, south latitude is negative, east longitude is positive, west longitude is negative.  
The variables refPtLat and refPtLong are the latitude and longitude of the reference point from which the distance and azimuth to the instantaneous pole of relative motion is to be determined.

refPtLat=30;

refPtLong=-130;

The variables rotPoleLat and rotPoleLong are the latitude and longitude of the pole around which the observed plate is moving in a positive (anticlockwise) direction.  These data can be obtained from the MORVEL plate-motion model of DeMets and others (2010).  For the motion of the Pacific Plate relative to the stable interior of the North American Plate, the pole is located at 48.9°S latitude (-48.9) and 108.3°E longitude (108.3).

rotPoleLat=80;

rotPoleLong=30;

The variable angVelDeg gives the angular velocity in degrees per million years.  For the motion of the Pacific Plate relative to the stable interior of the North American Plate, the angular velocity is 0.750°/Myr (Demets and others, 2010).

angVelDeg=0.750;

The variable earliestTime is the beginning of the time interval that is being modeled.  Past time is input as a negative number, future time is a positive number, and our unit of time is million years.  So an input value of -1 means 1 million years ago.

earliestTime=-4;

The variable latestTime is the end of the time interval being modeled.

latestTime=4;

The variable timeInterval defines the time steps that the model is to use in computing positions.

timeInterval=1;

How this model works mathematically

text

Some User - Defined Functions

The function makeVector computes the unit location vector of a given point on an assumed-spherical Earth, given that point's latitude and longitude.  The x coordinate is equal to Cos[latitude]*Cos[longitude].  The y coordinate is equal to Cos[latitude]*Sin[longitude].  The z coordinate is equal to Sin[latitude].  The input latitude and longitude values must be in decimal degrees (i.e., not degrees, minutes, seconds).  The forms "lat *(π/180)" and "long *(π/180)" in the square brackets convert the input data from degrees to radians, as required for the built-in Mathematica trigonometric functions (e.g., Sin, Cos, Tan, ArcSin, ArcCos, ArcTan).

makeVector[lat_,long_]:={Cos[lat*(π/180)]*Cos[long*(π/180)],Cos[lat*(π/180)]*Sin[long*(π/180)],Sin[lat*(π/180)]};

The function unitVector finds the corresponding vector with a length of 1 (i.e., unit length) given a vector of arbitrary length.  It uses a built-in Mathematica function Norm, which returns the length of the input vector.  Given a vector A with components {(Subscript[a, x]),  (Subscript[a, y]), (Subscript[a, z])}, the length or norm of A is equal to

Norm[A] = (Sqrt[((((Subscript[a, x])))^2)+((((Subscript[a, y])))^2)+((((Subscript[a, z])))^2)]).

The unit-vector that corresponds to A has components {((Subscript[a, x])/Norm[A]), ((Subscript[a, y])/Norm[A]), ((Subscript[a, z])/Norm[A])}.

unitVector[a_]:={a[[1]]/Norm[a],a[[2]]/Norm[a],a[[3]]/Norm[a]};

The function vectorAngle provides the angle between two vectors, in radians.  Recall that 360° =  (2 π) radians, so to convert from radians to degrees, multiply the radian value by (180/π).  
This function uses the vector dot product (also called an inner product or scalar product).   The notation "|A|" is the length of vector A, so if we note that vector A has components {(Subscript[a, x]),  (Subscript[a, y]), (Subscript[a, z])}, its length is

|A| ≡ Norm[A] = (Sqrt[((((Subscript[a, x])))^2)+((((Subscript[a, y])))^2)+((((Subscript[a, z])))^2)]).

Given two vectors of arbitrary length, A and B, the dot product of A and B is given by

A · B = |A| |B| Cos[θ] = ((Subscript[a, x]) (Subscript[b, x])) + ((Subscript[a, y]) (Subscript[b, y])) + ((Subscript[a, z]) (Subscript[b, z]) ).

Consequently,

θ = ArcCos[(A · B)/(|A| |B|)] = ArcCos[(((Subscript[a, x]) (Subscript[b, x])) + ((Subscript[a, y]) (Subscript[b, y])) + ((Subscript[a, z]) (Subscript[b, z]) ))/(|A| |B|)].

For more information, consult any textbook on vector analysis (e.g., Davis and Snider, 1987).

vectorAngle[a_,b_]:=ArcCos[Dot[a,b]/(Norm[a]*Norm[b])];

The functions vectorLat and vectorLong convert a Cartesian location vector back to a geographic coordinate system.

vectorLat[x_]:=ArcSin[x[[3]]]*(180/π);

vectorLong[x_]:=(x[[2]]/Abs[x[[2]]])*(ArcCos[x[[1]]/(Sqrt[(((x[[1]])^2))+(((x[[2]])^2))])])*(180/π);

The function angularDistance converts the location of two points, given in the geographic coordinate system of latitude and longitude, into Cartesian location vectors using the function makeVector, and then computes the angle between them using the function vectorAngle.

angularDistance[lat1_,long1_,lat2_,long2_]:=Module[{a,b,answer},a=makeVector[lat1,long1];b=makeVector[lat2,long2];answer=vectorAngle[a,b];answer];

The function zRotationMatrix enables computation of a rotation around the Z axis by an angle given by the product of the angular speed (w) and the change in time (Δt).

zRotationMatrix[w_,Δt_]:=(({{Cos[w Δt], -Sin[w Δt], 0}, {Sin[w Δt], Cos[w Δt], 0}, {0, 0, 1}}));

circularMotion[x_, w_, J1_, J2_, Δt_] := J1.zRotationMatrix[w, Δt].J2.x;

The function circMotionGeog takes a location vector (x), multiplies it by a coordinate transformation matrix (J2), a z - axis rotation matrix  (using function zRotationMatrix) through an angle given by the product of the angular speed (w) and the time interval (Δt), and the inverse coordinate transformation matrix (J1) , and converts the resulting location vector to geographic coordinates using the functions vectorLat and vectorLong.

circMotionGeog[x_, w_, J1_, J2_, Δt_, lat_, long_] := Module[{a, b, c, answer}, a = J1.zRotationMatrix[w, Δt].J2.x; b = vectorLat[a]; c = vectorLong[a]; answer = If[(Δt == 0), {lat, long, Δt}, {b, c, Δt}]; answer];

Define the matrices for coordinate transformation and rotation

First we define unit vectors along the right orthogonal Cartesian geographic coordinate axes, where the positive Z axis (zGeog) extends from Earth' s center toward the north pole, the X axis (xGeog) extends to the intersection of the Prime Meridian and the Equator, and the Y axis (yGeog) extends to longitude 90 °E on the Equator.

xGeog={1,0,0}; yGeog={0,1,0}; zGeog={0,0,1};

Second we define unit vectors along right orthogonal Cartesian axes in which the rotational axis around which the observed plate rotates in a positive (anti-clockwise) manner is the Z axis (zRot), the Y axis (yRot) coincides with the vector cross product zGeog × zRot = yRot, and the X axis (xRot) coincides with the vector cross product of yRot × zRot = xRot.

zRot=makeVector[rotPoleLat,rotPoleLong];

yRot=unitVector[Cross[zGeog,zRot]];

xRot=unitVector[Cross[yRot,zRot]];

Third we define the coordinate transformation matrix (transMatrix1) to transform a vector from the geographic Cartesian coordinate system (xGeog, yGeog, zGeog) to the rotational coordinate system (xRot, yRot, zRot).  This is a 3 x3 matrix of direction cosines between the different coordinate axes.

transMatrix1=(({{Cos[vectorAngle[xGeog,xRot]], Cos[vectorAngle[yGeog,xRot]], Cos[vectorAngle[zGeog,xRot]]}, {Cos[vectorAngle[xGeog,yRot]], Cos[vectorAngle[yGeog,yRot]], Cos[vectorAngle[zGeog,yRot]]}, {Cos[vectorAngle[xGeog,zRot]], Cos[vectorAngle[yGeog,zRot]], Cos[vectorAngle[zGeog,zRot]]}}));

Finally, we find the coordinate transformation matrix (transMatrix2) to transform a vector from the rotational coordinate system back to the geographic Cartesian coordinate system.

transMatrix2=Inverse[transMatrix1];

Adjust the input data

Convert the input angular speed (angVelDeg) from degrees per million years to radians per million years (ω)

ω=angVelDeg*(π/180);

Convert the geographic coordinates for the reference point (refPtLat, refPtLong) into a Cartesian unit location vector (refPtLocVector) using the function makeVector.

refPtLocVector=makeVector[refPtLat,refPtLong];

Output Data

The built-in Mathematica function Table is used to create a list of results, in which each record is composed of 3 data (latitude, longitude, and time)  enclosed in curly brackets.  The records span from earliestTime to latestTime in steps defined by the timeInterval, so the number of records is given by the result of ((latestTime - earliestTime) ÷ timeInterval) + 1.   The data in the table are generated by the function circMotionGeog, described above.

outputFile1=Table[circMotionGeog[refPtLocVector,ω,transMatrix2,transMatrix1,Δt,refPtLat,refPtLong],{Δt,earliestTime,latestTime,timeInterval}]

References Cited and Selected Bibliography

Cazenave, A., 1995, Geoid, topography and distribution of landforms, in Ahrens, T.J., [editor], Global Earth Physics -- A Handbook of Physical Constants:  American Geophysical Union Reference Shelf 1, p. 32-39, ISBN 0-87590-851-9.
Chase, C.G., 1972, The n-plate problem of plate tectonics:  Geophysical Journal of the Royal Astronomical Society, v. 29, p. 117–122.
Chase, C.G., 1978, Plate kinematics: the Americas, East Africa, and the rest of the world:  Earth and Planetary Science Letters, v. 37, p. 355–368.
Davis, H.F., and Snider, A.D., 1987, Introduction to Vector Analysis [5th edition]:  Boston, Massachusets, Allyn and Bacon, 365 p., ISBN 0-205-10263-8.
DeMets, C., Gordon, R.G., Argus, D.F. & Stein, S., 1990, Current plate motions:  Geophysical Journal International, v. 101, p. 425–478.
DeMets, C., Gordon, R.G., and Argus, D.F., 2010, Geologically current plate motions:  Geophysical Journal International, v.  121, p. x, doi: 10.1111/j.1365-246X.2009.04491.x
DeMets, C., Gordon, R.G., Argus, D.F., and Stein, S., 1994, Effect of recent revisions to the geomagnetic reversal timescale:  Geophysical Research Letters, v. 21, p. 2191–2194.
Minster, J.B., and Jordan, T.H., 1978, Present-day plate motions:  Journal of Geophysical Research, v. 83, p. 5331–5354.
Minster, J.B., Jordan, T.H., Molnar, P., and Haines, E., 1974, Numerical modeling of instantaneous plate tectonics:  Geophysical Journal of the Royal Astronomical Society, v. 36, p. 541–576.
Yoder, C.F., 1995, Astrometric and geodetic properties of Earth and the solar system, in Ahrens, T.J., [editor], Global Earth Physics -- A Handbook of Physical Constants:  American Geophysical Union Reference Shelf 1, p. 1-31, ISBN 0-87590-851-9.