Modeling-Afar

© 2019 by Vincent S. Cronin

Introduction

text

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

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

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

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=unitVect3D[Cross[z2,z1]];y2=unitVect3D[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.

The input parameters for the function tjCorner include...

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

cycloid[obsΩ_,movΩ_,refPt_,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=unitVect3D[obsΩ];e=Norm[obsΩ];f=unitVect3D[movΩ];g=Norm[movΩ];h=unitVect3D[d\[Cross]f];i=h;j=unitVect3D[h\[Cross]d];k=unitVect3D[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.refPt;answer];

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

distanceKM[cartRefPt_, cartOtherPt_] := Module[{radiusEarth, circumEarth, angDist, distance}, radiusEarth = 6371.02; circumEarth = 2*π*radiusEarth; angDist = VectorAngle[cartRefPt, cartOtherPt]*(180/π); distance = (angDist/360)*circumEarth; distance];

bearingAzimuth[cartRefPt_, cartOtherPt_] := Module[{vectorAlpha , northPole, vectorBeta, θ2, θ3, vectorGamma, θ4, bearing}, northPole = {0, 0, 1}; vectorAlpha = Cross[northPole, cartRefPt]; vectorBeta = Cross[cartOtherPt, cartRefPt]; θ2 = (180/π)*(VectorAngle[vectorAlpha, vectorBeta]); θ3 = (180/π)*(VectorAngle[vectorAlpha, cartOtherPt]); vectorGamma = Cross[refPtVect, vectorAlpha]; θ4 = (180/π)*(VectorAngle[vectorGamma, cartOtherPt]); bearing = If[((θ3 == 90) ∨ (θ3 == 180)), If[(θ4 <= 90), 0, 180], If[(θ3 > 90), (360 - θ2), θ2]]; bearing];

bearingDegDistanceKm[cartRefPt_, cartOtherPt_] := Module[{radiusEarth, circumEarth, angDist, distance, vectorAlpha , northPole, vectorBeta, θ2, θ3, vectorGamma, θ4, bearing, answer}, radiusEarth = 6371.02; circumEarth = 2*π*radiusEarth; angDist = VectorAngle[cartRefPt, cartOtherPt]*(180/π); distance = (angDist/360)*circumEarth; northPole = {0, 0, 1}; vectorAlpha = Cross[northPole, cartRefPt]; vectorBeta = Cross[cartOtherPt, cartRefPt]; θ2 = (180/π)*(VectorAngle[vectorAlpha, vectorBeta]); θ3 = (180/π)*(VectorAngle[vectorAlpha, cartOtherPt]); vectorGamma = Cross[refPtVect, vectorAlpha]; θ4 = (180/π)*(VectorAngle[vectorGamma, cartOtherPt]); bearing = If[((θ3 == 90) ∨ (θ3 == 180)), If[(θ4 <= 90), 0, 180], If[(θ3 > 90), (360 - θ2), θ2]]; answer = {"distance=",distance, "km, azimuth=",bearing, "degrees"}; answer];

13.3  Data from NNR-MORVEL56

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=(({{nb, 47.68, -68.44, 0.292}, {sm, 49.95, -84.52, 0.339}, {ar, 48.88, -8.49, 0.559}}))

hsWang=(({{nb, 40.2, 340.1, 0.2113}, {sm, 53.3, 318.3, 0.2248}, {ar, 37.3, 9.5, 0.5641}}));

nnrKreemer=(({{nb, 50.3991, -36.5919, 0.189645}, {sm, 54.6131, -71.5876, 0.225471}, {ar, 42.155, -0.270986, 0.456179}}));

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

eΩco = nnr[[1, 4]]*convert2Cart[nnr[[1, 2]], nnr[[1, 3]]];

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

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

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

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

eΩsa = nnr[[6, 4]]*convert2Cart[nnr[[6, 2]], nnr[[6, 3]]];

eΩan = nnr[[7, 4]]*convert2Cart[nnr[[7, 2]], nnr[[7, 3]]];

eΩjf = nnr[[8, 4]]*convert2Cart[nnr[[8, 2]], nnr[[8, 3]]];

eΩsm = nnr[[9, 4]]*convert2Cart[nnr[[9, 2]], nnr[[9, 3]]];

eΩar = nnr[[10, 4]]*convert2Cart[nnr[[10, 2]], nnr[[10, 3]]];

eΩin = nnr[[11, 4]]*convert2Cart[nnr[[11, 2]], nnr[[11, 3]]];

Argus

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

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

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

Kreemer

eΩnb = nnrKreemer[[1, 4]]*convert2Cart[nnrKreemer[[1, 2]], nnrKreemer[[1, 3]]];

eΩsm = nnrKreemer[[2, 4]]*convert2Cart[nnrKreemer[[2, 2]], nnrKreemer[[2, 3]]];

eΩar = nnrKreemer[[3, 4]]*convert2Cart[nnrKreemer[[3, 2]], nnrKreemer[[3, 3]]];

Wang

eΩnb = hsWang[[1, 4]]*convert2Cart[hsWang[[1, 2]], hsWang[[1, 3]]];

eΩsm = hsWang[[2, 4]]*convert2Cart[hsWang[[2, 2]], hsWang[[2, 3]]];

eΩar = hsWang[[3, 4]]*convert2Cart[hsWang[[3, 2]], hsWang[[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

text

caption

a=convert2Cart[13.209,40.091];aGeo=cart2geog[a];

b=convert2Cart[9.446,40.107];bGeo=cart2geog[b];

c=convert2Cart[11.317,41.442];cGeo=cart2geog[c];

d=convert2Cart[12.013,43.65];dGeo=cart2geog[d];

180-VectorAngle[Cross[a,c],Cross[c,b]]*(180/π)

180-VectorAngle[Cross[a,c],Cross[c,d]]*(180/π)

180-VectorAngle[Cross[c,b],Cross[d,c]]*(180/π)

inAR={a,c,d};

inNB={a,c,b};

inSM={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 3 million years of displacement.

modT=3;

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

outAR = Table[circMotion[inAR[[i]], eΩar, modT], {i, 1, Length[inAR]}];

displacedAR1=cart2geog[outAR[[1]]];

displacedAR2=cart2geog[outAR[[2]]];

displacedAR3=cart2geog[outAR[[3]]];

outNB = Table[circMotion[inNB[[i]], eΩnb, modT], {i, 1, Length[inNB]}];

displacedNB1=cart2geog[outNB[[1]]];

displacedNB2=cart2geog[outNB[[2]]];

displacedNB3=cart2geog[outNB[[3]]];

outSM = Table[circMotion[inSM[[i]], eΩsm, modT], {i, 1, Length[inSM]}];

displacedSM1=cart2geog[outSM[[1]]];

displacedSM2=cart2geog[outSM[[2]]];

displacedSM3=cart2geog[outSM[[3]]];

We use the function tjCorner to define the new position of the triple junction boundary for each plate (outTJX, outTJY and outTJZ).

outTJAR=tjCorner[outNB[[1]],outAR[[1]],outNB[[2]],outAR[[2]],outAR[[3]],outSM[[3]],outAR[[2]],outSM[[2]],outAR[[2]]]

displacedTJAR1=cart2geog[outTJAR[[1]]]

displacedTJAR2=cart2geog[outTJAR[[2]]]

displacedTJAR3=cart2geog[outTJAR[[3]]]

outTJSM = tjCorner[outNB[[3]], outSM[[1]], outNB[[2]], outSM[[2]], outAR[[3]], outSM[[3]], outAR[[2]], outSM[[2]], outAR[[2]]];

displacedTJSM1=cart2geog[outTJSM[[1]]]

displacedTJSM2=cart2geog[outTJSM[[2]]]

displacedTJSM3=cart2geog[outTJSM[[3]]]

outTJNB = tjCorner[outNB[[1]], outAR[[1]], outNB[[2]], outAR[[2]], outNB[[3]], outSM[[1]], outNB[[2]], outSM[[2]], outAR[[2]]];

displacedTJNB1=cart2geog[outTJNB[[1]]]

displacedTJNB2=cart2geog[outTJNB[[2]]]

displacedTJNB3=cart2geog[outTJNB[[3]]]

centroidDisplTJ=(outTJAR[[2]]+outTJSM[[2]]+outTJNB[[2]])/3;

cart2geog[centroidDisplTJ];

tjMigration=bearingDegDistanceKm[c,centroidDisplTJ]

VectorAngle[Cross[outTJAR[[1]],outTJAR[[2]]],Cross[outTJAR[[2]],outTJAR[[3]]]]*(180/π)

180-VectorAngle[Cross[outTJAR[[1]],outTJAR[[2]]],Cross[outTJAR[[2]],outTJAR[[3]]]]*(180/π)

VectorAngle[Cross[outTJSM[[1]],outTJSM[[2]]],Cross[outTJSM[[2]],outTJSM[[3]]]]*(180/π)

180-VectorAngle[Cross[outTJSM[[1]],outTJSM[[2]]],Cross[outTJSM[[2]],outTJSM[[3]]]]*(180/π)

VectorAngle[Cross[outTJNB[[1]],outTJNB[[2]]],Cross[outTJNB[[2]],outTJNB[[3]]]]*(180/π)

180-VectorAngle[Cross[outTJNB[[1]],outTJNB[[2]]],Cross[outTJNB[[2]],outTJNB[[3]]]]*(180/π)

We display the mid-point of each of the three polylines as a green dot

markers={outTJAR[[2]],outTJSM[[2]],outTJNB[[2]]};

out1=ListPointPlot3D[markers,AspectRatio->1,Boxed->False,PlotStyle->Directive[Green,PointSize[Medium]]];

and prepare the various polylines for graphic display.

out2= Graphics3D[{Brown,Line[outAR],Green,Line[outNB],Blue,Line[outSM],Red,Line[outTJAR],Line[outTJSM],Line[outTJNB],Black,Dashed,Line[initTJ]}, AspectRatio->1, Boxed->False];

Show[out2,out1]

caption.

VectorAngle[outTJAR[[2]],outTJSM[[2]]]*6371

VectorAngle[outTJAR[[2]],outTJNB[[2]]]*6371

VectorAngle[outTJSM[[2]],outTJNB[[2]]]*6371

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}, {"displacedTJAR1", displacedTJAR1[[1]], displacedTJAR1[[2]], modT}, {"displacedTJAR2", displacedTJAR2[[1]], displacedTJAR2[[2]], modT}, {"displacedTJAR3", displacedTJAR3[[1]], displacedTJAR3[[2]], modT}, {"displacedTJSM1", displacedTJSM1[[1]], displacedTJSM1[[2]], modT}, {"displacedTJSM2", displacedTJSM2[[1]], displacedTJSM2[[2]], modT}, {"displacedTJSM3", displacedTJSM3[[1]], displacedTJSM3[[2]], modT}, {"displacedTJNB1", displacedTJNB1[[1]], displacedTJNB1[[2]], modT}, {"displacedTJNB2", displacedTJNB2[[1]], displacedTJNB2[[2]], modT}, {"displacedTJNB3", displacedTJNB3[[1]], displacedTJNB3[[2]], modT}, {"displacedAR1", displacedAR1[[1]], displacedAR1[[2]], modT}, {"displacedAR12", displacedAR2[[1]], displacedAR2[[2]], modT}, {"displacedAR3", displacedAR3[[1]], displacedAR3[[2]], modT}, {"displacedNB1", displacedNB1[[1]], displacedNB1[[2]], modT}, {"displacedNB2", displacedNB2[[1]], displacedNB2[[2]], modT}, {"displacedNB3", displacedNB3[[1]], displacedNB3[[2]], modT}, {"displacedSM1", displacedSM1[[1]], displacedSM1[[2]], modT}, {"displacedSM2", displacedSM2[[1]], displacedSM2[[2]], modT}, {"displacedSM3", displacedSM3[[1]], displacedSM3[[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/Ryley/AfarWangHS.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.

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.

McKenzie, D.P., and Morgan, W.J., 1969, The evolution of triple junctions:  Nature, v. 224, p. 125-133.

Schouten, H., Cande, S.C., and Klitgord, K.D., 1985, Magnetic anomaly profiles, south (22°N to 28°N), in Rabinowitz, P.D., and Schouten, H., [editors], Mid-Atlantic ridge between 22° and 38°N:  Woods Hole, Massachusets, Marine Science International, Ocean Margin Drilling Program, Regional Atlas Series, Atlas 11.

Tucholke, B.E., and Schouten, H., 1988, Kane fracture zone:  Marine Geophysical Researches, v. 10, p. 1-39, 10 plates.

Vogt, P.R., 1986, Magnetic anomalies of the North Atlantic Ocean, in Vogt, P.R., and Tucholke, B.E., [editors], The western North Atlantic region:  Boulder, Colorado, Geological Society of America, The Geology of North America, volume M, plate 3.

Some web resources to investigate:

http://www.virtualocean.org/

http://www.geomapapp.org/

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