Triangle Strain in 2D

This Mathematica notebook was written by Vince Cronin in December, 2011.    

Revised January 3, 2012.

Introduction

The purpose of this notebook is to determine the horizontal strain in a triangular area between three non-colinear GPS stations, given their initial locations as well as their north-south and east-west velocities.  This code is based largely on the explanations published in Allmendinger, Cardozo and Fisher (2012) and Cardozo and Allmendinger (2009).

A stand-alone Macintosh application called SSPX performs this same analysis in a more comprehensive manner, including estimates of uncertainty.  SSPX is available for academic/research use for free via Nestor Cardozo’s website at http://homepage.mac.com/nfcd/work/programs.html.

Input data

GPS data from the EarthScope Plate Boundary Observatory is managed by UNAVCO (http://www.unavco.org/) and is available online for free at http://pbo.unavco.org/data.  The full public data holdings of UNAVCO are available via their “Data Archive Interface Version 2” at http://facility.unavco.org/data/dai2/app/dai2.html#.

For convenience, we compile the data in a 9 x 3 matrix in which each row corresponds to a different GPS station and the columns contain the following data types :  1 = site number, 2 = longitude, 3 = latitude, 4 = elevation (m, ellipsoid), 5 = UTM X coordinate, 6 = Y coordinate, 7  = N velocity (m/yr), 8 = E velocity, 9 = up velocity. 
The order in which the station data are presented in this matrix (i.e., which station corresponds to record 1, 2 or 3) is arbitrary.

inputData=(({{146, -120.537284, 39.337459, 2347.844, 712247.260, 4357118.219, 0.00715, -0.01025, -0.00181}, {149, -120.104975, 39.60212988, 2634.6887, 748566.261, 4387604.030, 0.00555, -0.00933, -0.00167}, {150, -120.033854, 39.29238024, 2619.0912, 755806.566, 4353418.468, 0.00626, -0.01044, -0.00357}}));

Initial computation

General explanation

We would like to use data from GPS arrays to define the average instantaneous strain rate in a triangular area between three GPS stations.  We know the initial location ((Subscript[X, 1]), (Subscript[X, 2])) of each GPS station from the station survey report or RINEX file.  We also know the east-west instantaneous velocity ((Subscript[u, 1])) and the north-south instantaneous velocity ((Subscript[u, 2])) of each station.  We do not know the elements of the deformation gradient tensor ((Subscript[e, 11]), (Subscript[e, 12]), (Subscript[e, 21]), (Subscript[e, 22])) or the coordinates of the translation vector ((Subscript[t, 1]), (Subscript[t, 2])).  The underlying relationships in matrix form is

(({{(^1)(Subscript[u, 1])}, {(^1)(Subscript[u, 2])}, {(^2)(Subscript[u, 1])}, {(^2)(Subscript[u, 2])}, {(^3)(Subscript[u, 1])}, {(^3)(Subscript[u, 2])}})) = (({{1, 0, (^1)(Subscript[X, 1]), (^1)(Subscript[X, 2]), 0, 0}, {0, 1, 0, 0, (^1)(Subscript[X, 1]), (^1)(Subscript[X, 2])}, {1, 0, (^2)(Subscript[X, 1]), (^2)(Subscript[X, 2]), 0, 0}, {0, 1, 0, 0, (^2)(Subscript[X, 1]), (^2)(Subscript[X, 2])}, {1, 0, (^3)(Subscript[X, 1]), (^3)(Subscript[X, 2]), 0, 0}, {0, 1, 0, 0, (^3)(Subscript[X, 1]), (^3)(Subscript[X, 2])}})) (({{(Subscript[t, 1])}, {(Subscript[t, 2])}, {(Subscript[e, 11])}, {(Subscript[e, 12])}, {(Subscript[e, 21])}, {(Subscript[e, 22])}}))

Following the explanation in Allmendinger and others (2012, p. 156), the matrix equation above has the form of

y = Mx

where y is a 1x6 matrix of known quantities (the instantaneous displacement/velocity vectors), x is a 1x6 matrix of unknown quantities (the translation vector and the deformation gradient tensor), and M is a 6x6 matrix of known values including zeros, ones and the location vector coordinates of the three GPS stations.  We need to rearrange this matrix equation so that all of the known quantities are collected on one side and the unknowns are on the other side of the equation.  This requires us to compute the inverse of M

x = (M^(-1))y

Mathematica has a built-in function to invert matrices, which we will employ in the form Inverse[M] =  (M^(-1)).  Allmendinger and others (2012) state that for perfectly constrained cases like this notebook was designed to handle, in which data from just three non-colinear GPS stations are used, the matrix can be inverted using LU decomposition.  Mathematica has a built-in process for LU decomposition.  For non-square matrices reflecting over-constrained cases with more GPS data than is minimally required to obtain a solution, Allmendinger and others (2012) recommend using the following formulation that employs transpose M matrices (after Press and others, 1986, and Menke, 1984)

x =[ (M^T)M] (M^T)y

Once the components of the deformation gradient tensor are known, we complete the 2-D process as follows.  Given the x matrix we just computed, we have the components of the displacement gradient tensor (Subscript[e, ij])

(Subscript[e, ij]) = (({{(Subscript[e, 11]), (Subscript[e, 12])}, {(Subscript[e, 21]), (Subscript[e, 22])}}))

The 2-D strain tensor is given by

(Subscript[ε, ij]) = (1/2)((Subscript[e, ij]) + (Subscript[e, ji]))  = (({{(Subscript[e, 11]), ((((Subscript[e, 12]) + (Subscript[e, 21])))/2)}, {((((Subscript[e, 21]) + (Subscript[e, 12])))/2), (Subscript[e, 22])}}))

and the 2-D rotation tensor is

(Subscript[ω, ij]) = (1/2)((Subscript[e, ij]) - (Subscript[e, ji])) = (({{0, ((((Subscript[e, 12]) - (Subscript[e, 21])))/2)}, {((((Subscript[e, 21]) - (Subscript[e, 12])))/2), 0}}))

The angle of rotation is equal to |(Subscript[ω, 12])| = |(Subscript[ω, 21])|.

Insert known data into governing equation to find unknowns

matrixM=(({{1, 0, inputData[[1,5]], inputData[[1,6]], 0, 0}, {0, 1, 0, 0, inputData[[1,5]], inputData[[1,6]]}, {1, 0, inputData[[2,5]], inputData[[2,6]], 0, 0}, {0, 1, 0, 0, inputData[[2,5]], inputData[[2,6]]}, {1, 0, inputData[[3,5]], inputData[[3,6]], 0, 0}, {0, 1, 0, 0, inputData[[3,5]], inputData[[3,6]]}}));

inverseM=Inverse[matrixM];

matrixY={{inputData[[1,8]]},{inputData[[1,7]]},{inputData[[2,8]]},{inputData[[2,7]]},{inputData[[3,8]]},{inputData[[3,7]]}};

In Mathematica, matrix multiplication is done using a dot symbol -- that is, a period -- or by using the built-in function Dot[m,k] where m and k are the names of the two matrices.

matrixX=Flatten[inverseM.matrixY];

The matrix eij is the 2-D displacement gradient tensor

eij={{matrixX[[3]],matrixX[[4]]},{matrixX[[5]],matrixX[[6]]}};

displGradTensor=MatrixForm[%];

The matrix epsij is the 2-D Lagrangian strain tensor

epsij={{eij[[1,1]],((eij[[1,2]]+eij[[2,1]])/2)},{((eij[[2,1]]+eij[[1,2]])/2),eij[[2,2]]}};

lagrangeStrainTensor=MatrixForm[%];

Rotation

omega12 is the rotation angle in degrees;  positive is a counterclockwise rotation.

omega12=((eij[[1,2]]-eij[[2,1]])/2)(180/π);

Translation

transCoord are the coordinates of the translation vector

transCoord={matrixX[[1]],matrixX[[2]]};

The length of the translation vector (transDistance) is given in meters.

transDistance=(Sqrt[((transCoord[[1]])^2)+((transCoord[[2]])^2)]);

unitTransVect is the unit vector that is coincident with the translation vector.

unitTransVect={(transCoord[[1]]/transDistance),(transCoord[[2]]/transDistance)};

transAngle is the angle between the north-directed vector and the translation vector.

transAngle=ArcCos[northUnitVector.unitTransVect](180/π);

The azimuth of the translation vector (transAzimuth) is measured clockwise from north.

transAzimuth=If[(unitTransVect[[1]]<0),(360-transAngle),transAngle];

Work with the greater and lesser principal axes

vectorNorm and unitVector2D are two user-defined functions that determine the length of a 2-D vector and find the unit vector corresponding to an arbitrary 2-D vector, respectively.

vectorNorm[x_]:=(Sqrt[x.x]);

unitVector2D[x_]:={x[[1]]/vectorNorm[x],x[[2]]/vectorNorm[x]};

northUnitVector is defined as a unit vector that points north, and from which azimuths can be computed.

northUnitVector={0,1};

eigenVects is the set of eigen vectors for the symmtrical Lagrangian strain tensor.  They are unit vectors that coincide with the principal axes of the strain tensor.

eigenVects=Eigenvectors[epsij];

eigenVector1={eigenVects[[1,1]],eigenVects[[1,2]]};

unitEVect1=unitVector2D[eigenVector1];

eigenVector2={eigenVects[[2,1]],eigenVects[[2,2]]};

unitEVect2=unitVector2D[eigenVector2];

The eigenvalues are the magnitudes of the principal strains

axisLengths=Eigenvalues[epsij];

The variable "a" is the length of the minimum extension axis of the 2-D horizontal strain ellipse;  "b" is the length of the maximum extension axis.

Variable a is the extension along the maximum extension axis, and variable a1 is the length of the corresponding semi - major axis after deformation of a circle that had an initial radius of 1.

a=If[(Abs[axisLengths[[1]]]>Abs[axisLengths[[2]]]),axisLengths[[1]],axisLengths[[2]]];

a1=1+a;

minExtAxis is the minimum extension axis, defined as being associated with the larger eigenvalue of the two axisLengths.

minExtAxis=If[(Abs[axisLengths[[1]]]>Abs[axisLengths[[2]]]),unitEVect1,unitEVect2];

angleA is the angle between the north-directed vector and the minimum extension axis.

angleA=ArcCos[northUnitVector.minExtAxis](180/π);

minExtAxisAz is the azimuth of the minimum extension axis.

minExtAxisAz=If[(minExtAxis[[1]]<0),(360-angleA),angleA];

Variable b is the extension along the maximum extension axis, and variable b1 is the length of the corresponding semi - major axis after deformation of a circle that had an initial radius of 1.

b=If[(Abs[axisLengths[[1]]]<Abs[axisLengths[[2]]]),axisLengths[[1]],axisLengths[[2]]];

b1=1+b;

maxExtAxis is the maximum extension axis, defined as being associated with the smaller eigenvalue of the two axisLengths.

maxExtAxis=If[(Abs[axisLengths[[1]]]>Abs[axisLengths[[2]]]),unitEVect2,unitEVect1];

angleB is the angle between the north-directed vector and the maximum extension axis.

angleB=ArcCos[northUnitVector.maxExtAxis](180/π);

maxExtAxisAz is the azimuth of the maximum extension axis.

maxExtAxisAz=If[(maxExtAxis[[1]]<0),(360-angleB),angleB];

Area strain

The area of a circle that has a radius of r is equal to π(r^2), and the area of an ellipse with semi-major and semi-minor axes a and b is equal to πab.

circleArea=π;

ellipseArea=π a1 b1;

areaStrain=(ellipseArea-circleArea)/circleArea;

Prepare to plot the results

The angle between the X axis (that is, the east-west axis) of the original GPS-station triangle and the longer principal axis is called "theta0".

theta0=If[(minExtAxis[[1]]<0),((π/180)(90-(minExtAxisAz-180))),((π/180)(90-minExtAxisAz))];

The angle theta combines theta0 with the rigid-body rotation omegaij so that the axes of the resulting plot image are north-south and east-west.

theta=theta0+(omega12*(π/180));

rotatedEllipse=Table[{(a1 Cos[t] Cos[theta])+(b1 Sin[t]*(-Sin[theta])),(a1 Cos[t] Sin[theta])+(b1 Sin[t] Cos[theta])},{t,0,2π,π/36}];

The plot range should be at least a bit larger than the strain ellipse, and this is handled by a variable called maxplot.

maxplot=1.5 a1;

The first plot file (plot1) plots a black ellipse with semimajor axis length of "a", semiminor axis length "b", and major axis inclined "theta" degrees relative to the east-directed vector.

plot1=ListLinePlot[rotatedEllipse, AspectRatio ->1, AxesLabel->{"East","North"},PlotRange->{{-maxplot,maxplot},{-maxplot,maxplot}}, PlotStyle->{Black}];

The second plot file (plot2) plots a dashed green reference circle.

referenceCircle=Table[{(Cos[t] Cos[theta])+(Sin[t]*(-Sin[theta])),(Cos[t] Sin[theta])+(Sin[t] Cos[theta])},{t,0,2π,π/36}];

plot2=ListLinePlot[%, AspectRatio ->1, PlotRange->{{-maxplot,maxplot},{-maxplot,maxplot}}, PlotStyle->{Thick,Dashed,Green}];

The third plot file (plot3) creates a blue line along the major axis of the ellipse.

plot3 = ListLinePlot[{{a1 Cos[theta], a1 Sin[theta]}, {-a1 Cos[theta], -a1 Sin[theta]}}, AspectRatio->1,PlotRange->{{-maxplot,maxplot},{-maxplot,maxplot}}, PlotStyle->{Blue}];

The fourth plot file (plot4) creates a red line along the minor axis of the ellipse.

plot4 = ListLinePlot[{{b1 Cos[theta + (π/2)], b1 Sin[theta + (π/2)]}, {-b1 Cos[theta + (π/2)], -b1 Sin[theta + (π/2)]}}, AspectRatio->1,PlotRange->{{-maxplot,maxplot},{-maxplot,maxplot}}, PlotStyle->{Red}];

The fifth plot file (plot5) merges the previous four plot files into a single file for display.

plot5=Show[plot2,plot1,plot3,plot4];

Results

Stations used

{inputData[[1,1]],inputData[[2,1]],inputData[[3,1]]}

Cartesian components of translation vector

transCoord

Length (meters) and Azimuth (degrees) of translation vector

{transDistance,transAzimuth}

Displacement gradient tensor

displGradTensor

Lagrangian strain tensor

lagrangeStrainTensor

Magnitude and azimuth of maximum principal extension

{b,maxExtAxisAz}

Magnitude and azimuth of maximum principal extension

{a,minExtAxisAz}

Area strain (dilation; negative strain indicates less area after deformation)

areaStrain

Infinitesimal rotation (clockwise + ve)

omega12

The plot of the strain ellipse (black) relative to the coordinate axes oriented north-south east-west (plot5) is below.  The axes are labeled in dimensionless strain.  The red axis is the minimum extension axis.  The green dashed circle shows a circle prior to deformation, so the black ellipse depicts that same circle after deformation.

plot5

References

Allmendinger, R.W., Loveless, J.P., Pritchard, M.E., and Meade, B., 2009, From decades to epochs -- Spanning the gap between geodesy and structural geology of active mountain belts:  Journal of Structural Geology, v. 31, p. 1409-1422, doi:  10.1016/j.jsg.2009.08.008.

Allmendinger, R.W., Cardozo, N, and Fisher, D.M., 2012, Structural Geology Algorithms, Vectors and Tensors :  Cambridge University Press, 289 p., ISBN 978 - 1 - 107 - 40138 - 9.

Cardozo, N., and Allmendinger, R.W., 2009, SSPX -- A program to compute strain from displacement/velocity data:  Computers and Geosciences, v. 35, p. 1343-1357, doi:  10.1016/j.cageo.2008.05.008.

Ferguson, John, 1994, Introduction to Linear Algebra in Geology :  London, Chapman & Hall, 203 p., ISBN 0 - 412 - 49350 - 0.

Fossen, Haakon, 2010, Structural Geology :  New York, Cambridge University Press, 463 p., ISBN 978 - 0 - 521 - 51664 - 8.

Menke, W., 1984, Geophysical  Data Analysis -- Discrete Inverse Theory:  Orlando, Florida, Academic Press.

Oertel, G., 1996, Stress and Deformation -- A Handbook on Tensors in Geology :  New York, Oxford University Press, 292 p., ISBN 0 - 19 - 509503 - 0.

Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes -- The Art of Scientific Computing:  Cambridge, Cambridge University Press.

Pollard, D.D., and Fletcher, R.C., 2005, Fundamentals of Structural Geology :  New York, Cambridge University Press, 500 p., ISBN 0 - 521 - 83927 - 0;  see especially Chapter 6.

Wolfram MathWorld has a web resource on eigenvectors and eigen decomposition available at 
http : // mathworld.wolfram.com/Eigenvector.html

Other Mathematica resources are available for download at 
Emmanuel Amiot, "Eigenvectors by Hand" from The Wolfram Demonstrations Project  http://demonstrations.wolfram.com/EigenvectorsByHand/
Yaroslav Bulatov, "Linear Transformation with Given Eigenvectors" from The Wolfram Demonstrations Project  http://demonstrations.wolfram.com/LinearTransformationWithGivenEigenvectors/