© 1988-2017 by Vincent S. Cronin;  revised March 2, 2017

6.5.1  Introduction

text.

6.5.2  Some user-defined functions

We will use some of the user-defined functions developed in previous chapters.

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

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

zRotation[w_,dT_]:={{Cos[(w dT)Degree],-Sin[(w dT)Degree],0},{Sin[(w dT)Degree],Cos[(w dT)Degree],0},{0,0,1}};

circMotion[x_,m1_,w_,dT_]:=Module[{m2,m3,answer},m2=zRotation[w,dT];m3=Inverse[m1];answer=m3.m2.m1.x;answer];

6.5.3  Instantaneous velocity vectors

Velocities are always measured relative to a specified frame of reference.  In this section, we are going to consider the velocity of a point P on one plate as observed from another plate.

Velocities are expressed in units of length per time;  angular velocities are in units of angle per time.  Strictly speaking, an instantaneous velocity is expressed in the limit as time shrinks to nearly zero.  The instantaneous velocity vector associated with a given point P is tangent to the surface of the sphere at that point.  Of course, there are an infinite number of vectors that are tangent to a sphere at a given point, so how do we determine the direction of the correct vector?

Let’s use the name pole for the unit location vector to the pole around which the plate containing point P moves in a positive direction (counter-clockwise) as observed from the other plate.  The current best source for pole locations and angular velocities for one plate as observed from another is the MORVEL-56 model of Chuck Demets and others (2010;  http://geoscience.wisc.edu/~chuck/MORVEL/).  The location vector to point P will be called ptP.

How can we find the direction of instantaneous motion at point P ?  Here are three ways to visualize the instantaneous velocity vector.  First, imagine a path that leads from point P directly to the pole, along a great-circle arc.  As you stand at P and look toward the pole, the instantaneous velocity vector extends to the right, perpendicular to that great-circle arc toward the pole.  Second, if you know the azimuth from a given point to the pole around which that point rotates in a positive direction, the direction of motion will be 90° greater than the azimuth to the pole.  Third, we can recognize that vectors pole and ptP are part of a unique plane through the center of the sphere, and the instantaneous velocity vector at point P is parallel to the vector result of pole × ptP, which is a vector normal to the pole - ptP plane.

What is the magnitude of the instantaneous tangential velocity vector at a given point?  How fast is a given point going?  All points on a rigid plate have the same angular velocity.  The magnitude of the tangential velocity of a given point (measured in units of length per time, such as km/Myr) is a function of the distance from the point to the pole.  The maximum magnitude occurs at an angular distance of 90° to the pole, and the magnitude is zero at the pole.  When our model Earth rotates by 1° on its axis, a point on the equator, 90° angular distance from the north pole, moves by (1/360) times the circumference in kilometers, or

maxDistPerDegree=((6371.01*2*π)/360)

kilometers.  The estimated mean radius of Earth (6371.01) limits the number of significant figures to two to the right of the decimal.  We can use a user-defined function round2DecRight  to perform that feat for us

round2DecRight[x_]:=N[Round[x*100]/100];

round2DecRight[111.195]

The user-defined function round2DecRight employs the built-in function Round, which follows the convention that a number of the form x.5 is rounded to x+1 when x is an odd number, and to x otherwise.  Consequently, 111.195 rounds to 111.20, and Mathematica (unhelpfully) drops the trailing zero.  That said, it is generally considered to be best practice to avoid rounding values in the middle of a string of computations, but rather to consider the appropriate number of significant figures at the end.

If we define θ as the angular distance to the rotational pole from our reference point, we can generalize the equation above so that it computes the distance traveled by any point on the model Earth surface during a rotation of 1°.

distPerDegree=(((6371.01*2*π)/360))Sin[θ];

The sine function produces a maximum value of 1 when θ = 90° and a value of 0 when θ = 0°, which is consistent with the manner in which the magnitude of the tangential velocity (the tangential speed) varies with angular distance from the pole.  Now, let’s add the final component, a variable angular rotation around the pole.  We will wrap-up the whole package into a user-defined function to find the magnitude of the tangential velocity vector (tangSpeed) and test it out with a rotation of 1° and θ = 90° to see if we get the correct answer that we computed above.

tangSpeed[rotAngle_,angDist_]:=(rotAngle*6371.01*2*π*Sin[angDist Degree])/360;

tangSpeed[1,90]

kilometers.  Good.  Our user-defined function gave us the expected angle for the one case that we have already computed.  Let’s see if our function works properly on a set of points ranging from pole (θ = 0°) to pole (θ = 180°) in steps of 30°.

Table[tangSpeed[1,i],{i,0,180,30}]

Peachy.  The results are symmetrical around the maximum value (111.195 km at θ = 90°) and drop to zero at both poles, with no unexpected sign errors.  Bueno.  So now we have some way of knowing the direction of relative motion at an arbitrary point, and the tangential speed of relative motion at that point.  But if we are standing at that point on a plate boundary, we might want our tangential speed to be expressed in something more meaningful to our human scale of perception -- cm/yr rather than km/Myr, for example.  A conversion is in order.  How many cm/yr is the equivalent of 1 km/Myr?

1 km = 1,000 m = 100,000 cm

1 Myr = 1,000,000 yr, so

1 km/Myr = 100,000 cm/1,000,000 yr = 0.1 cm/yr

For example, if we had an angular velocity of 1°/Myr and a point located 60° from the rotational pole, that point would have a tangential velocity of 96.30 km/Myr or (0.1*96.30) = 9.63 cm/yr, ignoring the significant figures.  (If your input data is good to 10 meters, your output data is not going to be good to 1/100 of a centimeter!  But I digress...)

Let’s do a worked example, and then practice with a homework exercise.

Worked example

One of the points along our imaginary little plate has an initial location of 26.6199°N latitude and 2.68982°E longitude.  The rotational pole is located at 65°N latitude and 15°W longitude, and the instantaneous angular velocity is 0.938°/Myr.

Question:  What is the azimuth (expressed in degrees) and magnitude of the tangential instantaneous velocity vector (expressed in cm/yr) at the specified reference point, and what are the coordinates of the unit vector that is parallel to that tangential velocity?

Computation and head scratching:  Back in chapter 2, we learned how to find the angular distance from one point to another, as well as the azimuth of one point as viewed from another.  First, let's convert the input geographic data into unit vectors.

pole=convert2Cart[65,-15];refPt=convert2Cart[26.6199,2.68982];

Second, we will find the angular distance between the two points, in degrees, using the built-in Mathematica functions VectorAngle and Degree, and call that angle ϕ

ϕ=VectorAngle[pole,refPt]/Degree;

Third, we will initialize the angular velocity ω and express it in degrees per million years.

ω=0.938;

Fourth, we find the tangential speed, employing the user-defined function tangSpeed

tangentialSpeed=tangSpeed[ω, ϕ];

in km/Myr, which corresponds to

tangSpeedCmYr=0.1*tangentialSpeed;

in cm/yr.  Fifth, we find the unit vector that is parallel to the direction of the tangential velocity vector.

unitVelVect=unitVect3D[Cross[pole,refPt]];

Now that we have done all of the easy stuff, we need to find the azimuth of the tangential velocity vector.  When you are standing at the reference point looking due north, the azimuth of the vector is found by measuring an angle on the imaginary horizontal plane that is tangent to Earth’s assumed-spherical surface at the very reference point you are standing on.  The angle is measured starting at north (0° azimuth) and continuing in a clockwise manner until you reach the direction of the velocity vector.  The azimuth will be a positive number less than 360°.  First, let’s tell your computer where the unit location vector pointing to the north pole is.

north={0,0,1};

Second, we find the normals to the planes through the reference point, the north pole, and Earth’s center (normal1) and through the  reference point, the unit velocity vector, and Earth’s center (normal2).

normal1=Cross[refPt,north];normal2=Cross[refPt,unitVelVect];

Third, we find the angle (γ, gamma) between the two normal vectors (normal1 and normal2)

γ=VectorAngle[normal1,normal2]/Degree;

Fourth, we find a unit vector (towardNorth) that is on the plane that the reference point shares with the north pole, so that when you are sitting on the reference point looking toward this new unit vector, you will be looking toward north.

towardNorth=unitVect3D[Cross[normal1,refPt]];

Fifth, we find the angle (ψ, psi) between the vectors towardNorth and normal2.

ψ=VectorAngle[towardNorth,normal2]/Degree;

Finally, we use ψ to help us determine whether the azimuth is γ or (360 - γ).

azimuth=If[(ψ>90),360-γ,γ];

Answers:  The azimuth of the tangential instantaneous velocity vector at the reference point is

azimuth

degrees.  The magnitude of the tangential instantaneous velocity vector at the reference point, expressed in centimeters per year, is

tangSpeedCmYr

The coordinates of the unit vector parallel to the tangential instantaneous velocity vector at the reference point are

unitVelVect

A little work now to make life easier later

Part 1.  A user-defined function for finding the azimuth

This business of finding an azimuth looks like a task we are likely to need to perform again, so let’s make it into a user-defined function findAzLocVect, so named because we are going to find the azimuth between two location vectors.

findAzLocVect[homePt_,otherPt_]:=Module[{a,b,c,d,e,f,answer},a={0,0,1};b=Cross[homePt,a];c=Cross[homePt,otherPt];d=VectorAngle[b,c]/Degree;e=unitVect3D[Cross[b,homePt]];f=VectorAngle[e,c]/Degree;answer=If[(f>90),360-d,d];answer];

Let’s see if this function gives us the same answer we obtained above.

findAzLocVect[refPt, unitVelVect]

Part 2.  A user-defined function for finding the tangential velocity vector and its magnitude

Now let’s make like a squirrel burying our acorns for a hungry day, and build a user-defined function that will take geographic data for the location of a reference point and a positive pole, along with the angular velocity around that pole expressed in degrees per million years, and tell us the orientation of the unit tangential-motion vector and the magnitude of the tangential velocity vector (the tangential speed) expressed in kilometers per million years.

firstTry[latPole_,longPole_,latPoint_,longPoint_,angVel_]:=Module[{a,b,c,d,e, answer},a=convert2Cart[latPole,longPole];b=convert2Cart[latPoint,longPoint];c=VectorAngle[a,b]/Degree;d=(angVel*6371.01*2*π*Sin[c Degree])/360;e=unitVect3D[Cross[a,b]];answer={d,e};answer];

And now we follow Ronald Reagan’s favorite Russian proverb, trust but verify.  We just computed the magnitude of the tangential vector in km/Myr

tangentialSpeed

in km/Myr, and the unit velocity vector

unitVelVect

Let’s see if we get the same results from firstTry

theAnswer=firstTry[65,-15,26.6199,2.68982,0.938]

Bueno, but it seems a bit cumbersome to have a mixed list of numbers as the output.  How do we extract data values from the list provided by the function we just wrote.  The list has two primary elements or dimensions.  If you want the value of the tangential speed in km/Myr, you would use

theAnswer[[1]]

because it is the first element of the list.  If you want the three components of the unit tangential-motion vector in a list, you would use

theAnswer[[2]]

because it is the second element of the list.  Finally, if you wanted the Z value of the vector, you would use

theAnswer[[2,3]]

because it is the third value in the second element of the list.  As a lazy man, it occurs to me that I can combine the tangential vector length with the unit vector coordinates to produce a vector that is the right length and pointed in the right direction.  So let’s modify our function to do this, and rename it in grand manner tangentialVelocityComputer

tangentialVelocityComputer[latPole_,longPole_,latPoint_,longPoint_,angVel_]:=Module[{a,b,c,d,e, answer},a=convert2Cart[latPole,longPole];b=convert2Cart[latPoint,longPoint];c=VectorAngle[a,b]/Degree;d=(angVel*6371.01*2*π*Sin[c Degree])/360;e=unitVect3D[Cross[a,b]]*d;answer=e;answer];

Let’s see if we get the same results from the tangentialVelocityComputer

anotherAnswer=tangentialVelocityComputer[65,-15,26.6199,2.68982,0.938]

Is this the same as the previous answer?  The length of the vector we just computed should be the same as the tangential speed we computed earlier.  We can just use the built-in Mathematica function Norm to find the length of that vector.

theAnswer[[1]]

Norm[anotherAnswer]

We can see that they are the same.  Excellent.  And the unit vector of the vector result of our tangentialVelocityComputer should be the same as the unit vector we computed earlier.

theAnswer[[2]]

unitVect3D[anotherAnswer]

The two results are the same.  Peachy.  We will store our new function tangentialVelocityComputer in a safe place till we need it.  But what if we already have the Cartesian coordinates of the rotational pole and our reference point?  Let’s modify our function some more to handle that situation, and call our new function tangVelCompFromCart.  Naw, that’s a really ugly name.  Let’s call it cupid, because Cupid worked with arrows.

cupid[pole_,point_,angVel_]:=Module[{a,b,c, answer},a=VectorAngle[pole,point]/Degree;b=(angVel*6371.01*2*π*Sin[a Degree])/360;c=unitVect3D[Cross[pole,point]]*b;answer=c;answer];

Let’s see if we get the same results from cupid

yetAnotherAnswer=cupid[pole,refPt,0.938]

Firme.  We will use cupid soon.

Boiling out the fat, leaving the meat

Let’s simplify all of the foregoing mess and list the code we need to solve a simple problem.  In fact, let’s re-solve the problem we just solved, so that we know our boiled-down code is correct.  And to be sure that we’re not cheating, we will begin by clearing all of the old answers.

ClearAll[pole,refPt,ω,tanVelVect,tanSpeed,azimuth,cupid,convert2cart,unitVect3D,findAzLocVect];

Question:  What is the azimuth (expressed in degrees) and magnitude of the tangential instantaneous velocity vector (expressed in km/Myr and cm/yr) at a reference point at 26.6199°N latitude 2.68982°E longitude given a rotational pole at 65°N latitude 15°W longitude and an angular velocity of 0.938°/Myr?  What are the coordinates of the tangential velocity vector?

Input

poleLat=65;poleLong=-15;angVel=0.938
refPtLat=26.6199;refPtLong=2.68982;

User-defined functions

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

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

cupid[pole_,point_,angVel_]:=Module[{a,b,c, answer},a=VectorAngle[pole,point]/Degree;b=(angVel*6371.01*2*π*Sin[a Degree])/360;c=unitVect3D[Cross[pole,point]]*b;answer=c;answer];

findAzLocVect[homePt_,otherPt_]:=Module[{a,b,c,d,e,f,answer},a={0,0,1};b=Cross[homePt,a];c=Cross[homePt,otherPt];d=VectorAngle[b,c]/Degree;e=unitVect3D[Cross[b,homePt]];f=VectorAngle[e,c]/Degree;answer=If[(f>90),360-d,d];answer];

Computation

pole=convert2Cart[poleLat,poleLong];refPt=convert2Cart[refPtLat,refPtLong];

tanVelVect=cupid[pole,refPt,angVel];

tanSpeed=Norm[tanVelVect];

azimuth=findAzLocVect[refPt,tanVelVect];

Answer

The magnitude of the tangential velocity vector is

tanSpeed

km/Myr, or

0.1*tanSpeed

cm/yr.  The azimuth of the tangential velocity vector is

azimuth

degrees measured clockwise from north.  The coordinates of the tangential velocity vector (km/Myr) are

tanVelVect

and, expressed in cm/yr are

0.1*tanVelVect

Exercise 6.5-1.  The little bridge along the Parkfield-Coalinga road near Parkfield, California, is one of the most famous bridges in the world, at least to seismologists.  Seismologists who get out of their offices and look at actual faults, that is.  The bridge has been rebuilt many times, in large part because the San Andreas fault passes between the northeast and southwest ends of the bridge.  According to Google Earth, the bridge is located at 35.895172°N latitude and 120.434657°W longitude.  The current best estimate for the location of the instantaneous pole of motion around which the Pacific plate rotates relative to North America is 48.9°S latitude and 108.3°E longitude, with an angular velocity of 0.750°/Myr (DeMets and others, 2010).
Write a Mathematica notebook that computes the direction (azimuth) and magnitude (cm/yr) of the instantaneous tangential motion of the Pacific plate relative to the North American plate at the Parkfield bridge.

Exercise 6.5-2.  UNAVCO maintains web-accessible datasets from geodetic GPS stations throughout the United States and elsewhere.  One of the GPS stations is called CARH (short for CARH_SCGN_CN2001), and is located at 35.88839°N latitude and 120.43082°W longitude, near the Parkfield bridge on the Pacific side of the San Andreas fault.  The motion of CARH relative to the North American reference frame (NAM08) is expressed as north-south, east-west, and up-down velocity vectors on the time-series graphs that are accessible at http://www.unavco.org/instrumentation/networks/status/pbo/overview/CARH.  The CARH data were accessed on March 2, 2017, and the NAM08 velocities were 26.51±0.19 mm/yr toward north, 20.21±0.17 mm/yr toward west, and 1.02±0.24 mm/yr up.  Using only the horizontal velocities, determine the direction (azimuth) and magnitude (cm/yr) of the motion of CARH relative to the stable cratonic interior of North America.  (You can do this with a calculator, the Pythagorean Theorem, and the SOH CAH TOA mnemonic to help solve right-triangle problems.)  How does the GPS velocity vector compare with your answer for Exercise 7-3?  Now do the same for PBO station HOGS, located SW of CARH at 35.86672°N latitude and 120.47950°W longitude, where the horizontal velocities were 29.30±0.06 mm/yr toward north, 21.91±0.08 mm/yr toward west, and 2.80±0.40 mm/yr down  (http://www.unavco.org/instrumentation/networks/status/pbo/overview/HOGS).  CARH is located quite near to the San Andreas fault trace, and HOGS is 5 km from the fault on the Pacific side.  Is there a difference in tangential (horizontal) velocities at the two PBO sites, and if so, form a hypothesis to account for this difference.  How might you test your hypothesis?

6.5.4  Plotting instantaneous tangential velocity vectors

This section is under construction as of 13 February 2012, and will be added as soon as the author becomes smart enough to know how to plot the velocity vector arrows correctly.

6.5.5  References and relevant texts

DeMets, C., Gordon, R.G., and Argus, D.F., 2010, Geologically current plate motions:  Geophysical Journal International, doi: 10.1111/j.1365-246X.2009.04491.x;  also see http://geoscience.wisc.edu/~chuck/MORVEL/