Modeling-Mendocino

© 2019-2021 by Vincent S. Cronin

Introduction

This code was used to generate Figure 7 from Cronin (2021), modeling the Mendocino triple junction in 1 million years.

13.2  User-defined functions

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

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

findGeogCoord[vect_]:=Module[{lat,long,a,b,θ},a=ArcSin[vect[[3]]];θ=ArcCos[vect[[1]]/(Sqrt[(((vect[[1]])^2)+((vect[[2]])^2))])];b=If[(vect[[2]]<0),(-θ),(θ)];lat=a(180/π);long=If[((Abs[vect[[1]]]<(1×(10^(-14))))&&(Abs[vect[[2]]]<(1×(10^(-14))))),0,(b(180/π))];{lat,long}];

makeGreatCircle[normal_]:=Module[{a,x1,y1,z1,x2,y2,z2,j1,j2,j3,θ,b},a=Table[{Cos[i Degree],Sin[i Degree],0},{i,0,360,5}];x1={1,0,0};y1={0,1,0};z1={0,0,1};z2=normal;x2=Normalize[Cross[z2,z1]];y2=Normalize[Cross[z2,x2]];j1={{x1.x2,y1.x2,z1.x2},{x1.y2,y1.y2,z1.y2},{x1.z2,y1.z2,z1.z2}};θ=VectorAngle[z1,z2];j2={{1,0,0},{0,Cos[θ],Sin[θ]},{0,-Sin[θ],Cos[θ]}};j3=Inverse[j1];b=Table[j3.j2.j1.a[[i]],{i,1,Length[a]}];b]

Constructing the new ridge-ridge-ridge triple junction takes some work, so we created a function called tjCorner to help.  The problem that we confront with a RRR triple junction that we have not confronted before is that the ridge segments must propagate or grow during displacement.  There are a couple of ways of handling this issue;  here, we will maintain the ridge axes as great-circle arcs and create new 2-segment polylines to represent the new boundaries of each of the three plates.

tjCorner[seg1enda_,seg1endb_,seg1mida_,seg1midb_,seg2enda_,seg2endb_,seg2mida_,seg2midb_,tj0_]:=Module[{a,b,c,d,e,f,g,h,i,answer},a=Normalize[(seg1enda+seg1endb)/2];b=Normalize[(seg1mida+seg1midb)/2];c=Normalize[(seg2mida+seg2midb)/2];d=Normalize[(seg2enda+seg2endb)/2];e=a\[Cross]b;f=c\[Cross]d;g=Normalize[e\[Cross]f];h=VectorAngle[g,tj0];i=If[(h<(π/2)),g,-g];answer={a,i,d};answer];

cycloid[refΩ_,obsΩ_,obsPt_,t_]:=Module[{a,b,c,d,e,f,g,h,i,j,k,l,m1, m2, m3, m4, m5,n,p,answer},a={1,0,0};b={0,1,0};c={0,0,1};d=Normalize[refΩ];e=Norm[refΩ];f=Normalize[obsΩ];g=Norm[obsΩ];h=Normalize[d\[Cross]f];i=h;j=Normalize[h\[Cross]d];k=Normalize[i\[Cross]f];m1={{k[[1]],k[[2]],k[[3]]},{i[[1]],i[[2]],i[[3]]}, {f[[1]],f[[2]],f[[3]]}};m2={{Cos[(-g Degree) t],Sin[(-g Degree) t],0},{-Sin[(-g Degree) t],Cos[(-g Degree) t],0},{0,0,1}};l=VectorAngle[d,k];If[l<(π/2),(n=1),(n=(-1))];p = (n) (VectorAngle[f,d]);m3={{Cos[p],0,-Sin[p]},{0,1,0},{Sin[p],0,Cos[p]}};m4={{Cos[(e Degree) t],Sin[(e Degree) t],0},{-Sin[(e Degree) t],Cos[(e Degree) t],0},{0,0,1}};m5={{j[[1]],h[[1]],d[[1]]},{j[[2]],h[[2]],d[[2]]},{j[[3]],h[[3]],d[[3]]}};answer=m5.m4.m3.m2.m1.obsPt;answer];

circMotion[x_,angVelVect_,dT_]:=Module[{north,rotPole,w,xPole,yPole,m1,m2,m3,answer},north={0,0,1};rotPole=Normalize[angVelVect];w=Norm[angVelVect];xPole=Normalize[rotPole\[Cross]north]; yPole=Normalize[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 module findDistAzimuthVectInput[]) acts on unit location vectors to determine the distance (in km) and azimuth (in degrees) from a reference point (location vector refPtVect) to a target point (location vector otherPtVect).

findDistAzimuthVectInput[refPtVect_,otherPtVect_]:=Module[{radiusEarth, circumEarth, angDist, distance, northPole,vectorA,vectorB,θ1,θ2,vectorC,θ3,azimuth,answer},radiusEarth = 6371.02; circumEarth = 2*π*radiusEarth; angDist = VectorAngle[refPtVect,otherPtVect]*(180/π);distance = (angDist/360)*circumEarth; northPole={0,0,1};vectorA = Cross[northPole,refPtVect];vectorB = Cross[otherPtVect,refPtVect];θ1=VectorAngle[vectorA,vectorB](180/π);θ2=VectorAngle[vectorA,otherPtVect](180/π);vectorC=Cross[refPtVect,vectorA];θ3=VectorAngle[vectorC,otherPtVect](180/π);azimuth=If[((θ2==90)∨(θ2==180)),If[(θ3<=90),0,180],If[(θ2>90),(360-θ1),θ1]];answer = {"distance=",distance, "km, azimuth=",azimuth, "degrees"}; answer];

13.3  Data from external reference frame

As in chapter 12, we will use data from the no-net-rotation solution derived by Argus and others (2011) based on the MORVEL velocity model (DeMets and others, 2010) and Peter Bird’s descriptions of 56 plates (Bird, 2003).  The abbreviations in the first column are an=Antactica, co=Cocos, nz=Nazca, na=North America, nb=Nubia (western Africa), pa=Pacific, and sa=South America.  For a given row/record, the second column contains the latitude of the poles around which the corresponding plate moves in a counter-clockwise manner, the third column has the corresponding polar longitude, and the fourth column has the angular speed in degrees per million years.

nnrArchive=(({{co, 26.93, -124.31, 1.198}, {nz, 46.23, -101.06, 0.696}, {na, -4.85, -80.64, 0.209}, {nb, 47.68, -68.44, 0.292}, {pa, -63.58, 114.70, 0.651}, {sa, -22.62, -112.83, 0.109}, {an, 65.42, -118.11, 0.250}, {jf, -38.31, 60.04, 0.951}, {sm, 49.95, -84.52, 0.339}, {ar, 48.88, -8.49, 0.559}, {in, 50.37, -3.29, 0.544}}));

nnrArgus=(({{jf, -38.31, 60.04, 0.951}, {pa, -63.58, 114.70, 0.651}, {na, -4.85, -80.64, 0.209}}));

hsWang=(({{jf, -36.7, 59.7, 1.1191}, {pa, -59.3, 96.1, 0.7703}, {na, -35.3, 325.7, 0.1679}}));

We use these data to compute the angular velocity vector for each plate relative to the no-net-rotation external reference frame.

eΩjf=nnrArgus[[1,4]]*convert2Cart[nnrArgus[[1,2]],nnrArgus[[1,3]]];

eΩpa=nnrArgus[[2,4]]*convert2Cart[nnrArgus[[2,2]],nnrArgus[[2,3]]];

eΩna=nnrArgus[[3,4]]*convert2Cart[nnrArgus[[3,2]],nnrArgus[[3,3]]];

13.5  Modeling a triple junction

The triple junctions imagined in first-generation plate tectonics were simple geometric forms.  Many, most or all actual triple junctions are interestingly complex.  We will look at real triple junctions in a future chapter.  Here, we simply indicate how we can model a couple of triple junctions.

Ridge-ridge-ridge triple junction

Variable a is associated with the point on the trench between Juan de Fuca and North America

a=convert2Cart[42.285,-125.290];aGeo=cart2geog[a];

Variable b is associated with the point on the transform between Juan de Fuca and Pacific

b=convert2Cart[40.422,-126.869];bGeo=cart2geog[b];

Variable c is associated with the triple junction

c=convert2Cart[40.311,-124.416];cGeo=cart2geog[c];

Variable d is associated with the point on the transform between Pacific and North America

d=convert2Cart[38.541,-123.305];dGeo=cart2geog[d];

inNA={a,c,d};

inJF={a,c,b};

inPA={b,c,d};

We create a list of location vectors along the triple junction prior to displacement (initTJ).

initTJ={a,c,d,c,b};

We choose to examine the triple junction after a model time (modT) of 1 million years of displacement.

modT=1;

Rotate the boundary polylines for (modT) million years around its axis (eΩx) to its new position (outX).

outNA = Table[circMotion[inNA[[i]], eΩna, modT], {i, 1, Length[inNA]}];

displacedNA1=cart2geog[outNA[[1]]];

displacedNA2=cart2geog[outNA[[2]]];

displacedNA3=cart2geog[outNA[[3]]];

outJF = Table[circMotion[inJF[[i]], eΩjf, modT], {i, 1, Length[inJF]}];

displacedJF1=cart2geog[outJF[[1]]];

displacedJF2=cart2geog[outJF[[2]]];

displacedJF3=cart2geog[outJF[[3]]];

outPA = Table[circMotion[inPA[[i]], eΩpa, modT], {i, 1, Length[inPA]}];

displacedPA1=cart2geog[outPA[[1]]];

displacedPA2=cart2geog[outPA[[2]]];

displacedPA3=cart2geog[outPA[[3]]];

Finding the TJ

v1=Normalize[(outPA[[1]]+outJF[[3]])/2]

v2=Normalize[(outPA[[2]]+outJF[[2]])/2]

v3=Normalize[(outPA[[3]]+outNA[[3]])/2]

v4=Normalize[(outPA[[2]]+outNA[[2]])/2]

normalJFPA=Cross[Cross[v1,v2],Cross[v3,v4]];

newTJ = If[(VectorAngle[normalJFPA,outPA[[2]]]<(π/4)),normalJFPA,(-1*normalJFPA)];

answer=findDistAzimuthVectInput[c,newTJ];

N[%]

out2= Graphics3D[{Brown,Line[outNA],Green,Line[outJF],Blue,Line[outPA],Black,Dashed,Line[initTJ]}, AspectRatio->1, Boxed->False];

Show[out2]

caption.

Perhaps this mismatch is part of an explanation for the observed development of a tiny plate in the middle of this triple junction:  the Galapagos plate.

Make a table of your results for output to an Excel spreadsheet.  The first record contains headers. The next four records will be the four points that define the current triple junction.  Records 6 through 14 will be the nine points that define the displaced triple junction after modT million years.  Records 15 through 17 will be the three points that define the edge of the Arabian Plate displaced at modT million years.  Records 18 through 20 will be the three points that define the edge of the Somalian Plate displaced at modT million years.  Records 21 through 23 will be the three points that define the edge of the Nubian Plate displaced at modT million years.

resultsTable=(({{"plate", "latitude", "longitude", "MyrBP"}, {"TJ_boundary_pt", aGeo[[1]], aGeo[[2]], 0}, {"TJ_boundary_pt", bGeo[[1]], bGeo[[2]], 0}, {"TJ_boundary_pt", cGeo[[1]], cGeo[[2]], 0}, {"TJ_boundary_pt", dGeo[[1]], dGeo[[2]], 0}, {"displacedNA1", displacedNA1[[1]], displacedNA1[[2]], modT}, {"displacedNA2", displacedNA2[[1]], displacedNA2[[2]], modT}, {"displacedNA3", displacedNA3[[1]], displacedNA3[[2]], modT}, {"displacedJF1", displacedJF1[[1]], displacedJF1[[2]], modT}, {"displacedJF2", displacedJF2[[1]], displacedJF2[[2]], modT}, {"displacedJF3", displacedJF3[[1]], displacedJF3[[2]], modT}, {"displacedPA1", displacedPA1[[1]], displacedPA1[[2]], modT}, {"displacedPA2", displacedPA2[[1]], displacedPA2[[2]], modT}, {"displacedPA3", displacedPA3[[1]], displacedPA3[[2]], modT}}));

MatrixForm[resultsTable]

resultsTable

Export data

1.  In the same folder or directory in which this notebook is located, create a blank Excel spreadsheet and call it “Rutabega.xlsx”

2.  Enter the following as an input line (under the black box, below):

Export[,outFile1];

3. Put your cursor between the left square bracket and the comma in the Export[,outFile1] statement and click to establish the insertion point we will need in the next step.

4. Go to the Insert menu, select File Path, navigate to the correct input data file and choose it, and the correct file path will be inserted at the cursor.  In this example, the Excel file is located on the desktop of Vince’s computer, and so Mathematica will insert the path "/Users/vince/Desktop/Rutabega.xlsx" between the brackets.  The path will be different for every different file, and on every different computer.

Export[“/Users/vince/Desktop/Rutabega.xlsx”,];

You can make sure it is an input line by clicking on that line’s bracket on the far right edge of this window, going to the Format menu, selecting the Style submenu, and then choosing Input

The file created in the next line is an Excel (xlsx) file that will be placed in the same directory on your computer as this notebook file.

IMPORTANT NOTE:  The user must supply a name for the output file that is different from existing file names, or else the existing files may be over-written.

Export["/Users/vince/Desktop/Mendocino.Argus.xlsx", resultsTable];

13.7  References and some relevant texts

Argus, D.F., Gordon, R.G., and DeMets, C., 2010, 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., 1987, The contribution of cycloid wobble to plate boundary strain: Geological Society of America, Abstracts with Programs, v. 19, p. 149.

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., 1992, Types and kinematic stability of triple junctions: Tectonophysics, v. 207, p. 287-301.

Cronin, V.S., 1994, Instantaneous velocity of mid-ocean ridges: Tectonophysics, v. 230, p. 151-159.

Cronin, V.S., 2021, Triple junctions, in Alderton, D., and Elias, S.A., Encyclopedia of Geology [2nd edition]:  Elsevier, p. 947-956.

Klein, E.M., Smith, D.K., Williams, C.M., and Schouten, H., 2005, Counter-rotating microplates at the Galapagos triple junction:  Nature, v. 433, 24 February, p. 855-858.