Скачать 89.5 Kb.

Primer on infinitesimal strain analysis in 1, 2 and 3D The purpose of this document is to provide some background for the infinitesimal strain analysis we are about to undertake, making use of velocity data from GPS sites. It is best to study this document alongside the treatment of strain in an introductory structural geology textbook. 1dimensional analysis Imagine a thin strip of ideal elastic material – an ideal rubber band, if you will. You measure its initial length, and use the symbol l_{original} to represent that length. Then you stretch it, measure the stretched strip, and call its new length l_{final}. Two descriptive terms are particularly useful to us in describing your observations. The elongation or extension of the strip (e) is defined as extension = e = If the strip is longer after you change its length, the extension is a positive number, and it is negative if you shorten the strip. If there is no change in length, the extension is zero. Extension is a dimensionless number, because it is the quotient of a length divided by a length. The stretch (S) of the strip is defined as stretch = S = e + 1 = l_{final}/l_{original} A positive extension results in a stretch that is greater than 1, and a negative extension or shortening results in a stretch that is less than 1. Extension (or elongation) and stretch are terms that an introductory student of structural geology needs to fully comprehend. Now let’s imagine doing a similar experiment, in which we have an ideal, thin elastic strip whose end is fixed to the “0” on one end of a ruler. We will also imagine that the ruler is along the x axis of a Cartesian coordinate system whose origin is the “0” on the ruler. Before we change the length of the strip, we put a black mark on the elastic strip at a distance of 5 cm from the fixed end, and a red mark at a distance of 8 cm from the end (Figure 1). Figure 1. Change in length of an ideal elastic strip, showing the strip’s original and final lengths. The x coordinate axis extends parallel with the ruler. Subscripts: b = black, r = red, o = original, f = final The black mark has coordinates in its original position of {5, 0, 0}, but because this is a 1dimensional problem we can simplify and use just the x coordinate (x_{bo}): {5}. In a similar manner, the red mark has an original location (x_{ro}) of {8}. Now, we stretch the band so that the black mark is at 6 cm along the ruler, so its final position (x_{bf}) is {6}. The final length divided by the initial length is the stretch (S) – the very same stretch that we just defined above. S = l_{final}/l_{original} The extension (call it e_{b}) and stretch (call it S_{b}) associated with the black mark is given by e_{b} = (6 – 5)/5 = 0.2 S_{b} = e_{b} + 1 = 1.2 or S_{b} = 6/5 = 1.2 The stretch will be the same for all points along the elastic strip, including the black and red marks. The stretch is also known as the 1dimensional Lagrangian deformation gradient in the x coordinate direction. Without digressing into history or biography, a Lagrangian approach to these problems (named after JosephLouis Lagrange) starts with an initial configuration and projects to a final, deformed configuration, while an Eulerian approach (named after Leonhard Euler, whose last name sounds like “oiler”) starts with the deformed and projects backward to the original. We are only going to use the Lagrangian approach in the GPS strain problem. If we know the initial coordinate {x_{original}} of an arbitrary point along the thin elastic strip, we can compute the final coordinate of that point {x_{final}} given the deformation gradient or stretch (S). Elongating the strip so that the black mark moves from 5 to 6 on the ruler will result in a deformation gradient or stretch of (6/5) = 1.2. That same elongation will cause the position of the red mark to change from 8 to (1.2 x 8) = 9.6 if the ideal strip we are stretching is homogeneous and has a linearelastic rheology. Another way of defining the deformation gradient focuses simply on the change in distance between the two marks on our elastic strip. The initial distance between the black and red marks (x_{o}) is 8 – 5 = 3, while the final distance between the black and red marks (x_{f}) is 9.6 – 6 = 3.6. The stretch between these two points (S) is S = The stretch is the slope of a line that relates the difference in initial position of two points to the difference in position of the same two points at some other time. In the case involving our black and red marks, S = 3.6/3 = 1.2 and so the deformation gradient is 1.2. For our purposes as we build toward the GPS strain problem, it is useful to think about the 1D changes in the elastic strip in terms of displacements, which are commonly referenced using the symbol u. Let us define the location vector to the initial location of the black mark as x_{bo}, whose 1D coordinate is {5}, and the location vector to the final location of the black mark is x_{bf}_{ }= {6}. The displacement vector u_{b} is given by the vector equation u_{b} = x_{bf} – x_{bo}, so u_{b} = {6} – {5} = {1} This is called a displacement vector, because it extends from the original position to the displaced or final position. Similarly, let’s define the location vector to the initial location of the red mark as x_{ro}, whose 1D coordinate is {8}, and the location vector to the final location of the red mark is x_{rf}_{ }= {9.6}. The displacement vector u_{r} is given by the vector equation u_{r} = x_{rf} – x_{ro}, so u_{r} = {9.6} – {8} = {1.6} The displacement vector associated with the black mark that was originally located at {5} has a length of 1, while the displacement vector associated with the red mark initially 3 cm farther down the ruler at {8} has a length of 1.6. The magnitude of the displacement vector varies as a function of position along the x coordinate axis. Earlier, we defined the extension (e) as extension = e = The extension is also equal to the difference between the lengths of two displacement vectors divided by the difference between the lengths of the corresponding two initiallocation vectors e = . The extension is the slope of a line that relates the difference in initial position of two points to the difference in displacement of the same two points at some other time, or the 1dimensional displacement gradient. Using the 1D example involving the black and red marks along the elastic strip, The 1D deformation gradient is equal to the stretch, and in this example has a value of 1.2. The 1D displacement gradient is equal to the extension, and in this example has a value of 0.2. So far, we have thought about deformation and displacement that we can see and easily measure. The change in length of the strip between the origin and the black mark amounted to 20% of the original length. This kind of change is normally considered finite strain. The obvious strain that we observe in rocks is the result of finite strain. If the change is on the order of 1% or less of the original length, it would normally be considered to be infinitesimal strain. For example, a change of 5 millimeters in the distance between two GPS sites that were originally 35 km apart is an excellent example of infinitesimal strain, because the change would be ~1.4x10^{5 }%. In many cases this strain is ongoing, so that it makes more sense to consider a displacement rate or velocity rather than a single displacement. When considering velocity, which has dimensions of length per time, the change in length becomes infinitesimal as the corresponding time interval is reduced toward zero. Thus when time is involved, the idea of instantaneous strain or instantaneous velocity is comparable to infinitesimal strain. The infinitesimal 1D deformation gradient (stretch) can be expressed as and the infinitesimal 1D displacement gradient (extension) is The symbol indicates a partial derivative, which is not really necessary in the 1D case but prepares us to think about the 2D and 3D cases in which the functions might vary in 2 or 3 dimensions. 3dimensional analysis The material that follows is generally confined to information about infinitesimal strain analysis that is applicable to the analysis of GPS velocity data. More comprehensive treatments of this material are available from several sources, including Malvern (1969), Timoshenko and Goodier (1970), Means (1976), Oertel (1996), Pollard and Fletcher (2005), and Allmendinger and others (2012). If the mental image we worked with in the previous section was that of an ideal elastic string whose length is changed in one dimension, parallel to the x coordinate axis, we are now going to build on that experience to think of a continuous elastic sheet that can be deformed in 3 spatial dimensions with Cartesian x, y and z axes. Analogous to the 1dimensional case, the 3D Lagrangian displacement equation includes the displacement vector {u_{x}, u_{y}, u_{z}} associated with a particular reference point, a vector defining the initial location of the reference point {x_{o}, y_{o}, z_{o}}, and the 3x3 displacement gradient tensor. The matrix of extensions in the displacement gradient tensor is . The elements of the matrix along the diagonal of the displacement gradient tensor (e_{11}, e_{22}, e_{33}) are the extensions along the x, y and z coordinate axes, respectively. In other words, the diagonal terms are related to changes in length along each of three orthogonal directions. The offdiagonal or offaxis elements are the tensor shear strains of lines originally parallel to the coordinate axes, and reflect angular changes within the deforming material. (The tensor shear strain is half of the engineering shear strain you might have encountered elsewhere.) Component e_{21} expresses the counterclockwise (positive) rotation of a vector that is initially parallel to the x coordinate axis (or 1 axis) toward the y axis (or 2 axis) around the z axis (or 3 axis) during deformation. Similarly, given a vector that is initially parallel to the y axis, component e_{12} expresses the rotation of that vector clockwise toward the x axis around the z axis during deformation. The displacement gradient tensor is an asymmetric matrix that represents both distortion (that is, change in shape or change in volume) and rotation. A tensor that represents only rotation would be antisymmetric, so that e_{12} = e_{21}, e_{13} = e_{31}, and e_{23} = e_{32}. The asymmetry in the displacement gradient tensor arises because of infinitesimal rotations and distortion. An asymmetric tensor can be decomposed into the sum of a symmetric tensor and an antisymmetric tensor. The asymmetric displacement tensor e_{ij} shown above is the sum of the symmetric infinitesimal strain tensor _{ij} and the antisymmetric rotation tensor _{ij} : e_{ij} = _{ij} + _{ij}, where The antisymmetric rotation tensor _{ij} is also known as an axial vector. The axial vector is directed along the rotational axis, and the length of the axial vector is equal to the angle of rotation expressed in radian measure. The Cartesian coordinates of the axial vector are given by {r_{x}, r_{y}, r_{z}}, where _{.} The angle of rotation in radians is the length of the axial vector, If the angle is a positive number, the rotation is an anticlockwise rotation as viewed looking down the rotational axis toward the origin of the coordinate system. The computed angle only has meaning under conditions of plane strain, and reflects the angle and direction that the maximum principal strain axis is rotating as deformation progresses. If deformation is entirely by pure strain (also known as irrotational strain or coaxial strain), the angle will be zero. Example. Given the following displacement gradient tensor (e_{ij}), calculate the strain tensor (_{ij}) and the rotation tensor (_{ij}). Solution. The strain tensor is and the rotation tensor is The Cartesian coordinates of the axial vector are or {1, 2, 3}. The angle of rotation is radians Eigenvectors, eigenvalues, principle strains and shear strains The principal strain axes are found by computing the eigenvectors of the strain tensor _{ij}. The eigenvectors are vectors in the directions of the principal strain axes. The eigenvalues of _{ij} are the principal extensions in the principal directions. The largest of the three eigenvalues is the greatest principal extension, e_{1}, and the length of the semimajor axis of the strain ellipsoid is equal to the greatest principal stretch, S_{1} = e_{1} + 1. The semiminor axis is S_{3} = e_{3} + 1, where e_{3} is the smallest eigenvalue and the least principal extension, and the intermediate axis is S_{2} = e_{2} + 1, where e_{2} is the intermediate principal extension. The magnitude of the infinitesimal engineering shear strain (_{max}) at 45° to the maximum principal strain axis can be computed several ways. _{max} = _{max} = e_{1} – e_{3} and, recalling that (S_{1})^{2} = (e_{1} +1)^{2} and (S_{3})^{2} = (e_{3} +1)^{2}, _{max} =. The engineering shear strain really only has meaning in cases of plane strain. The volume strain, which is also the first invariant of the strain tensor, is given by volume strain = _{11} + _{22} + _{33} = e_{1} + e_{2} + e_{3} = (S_{1} S_{2} S_{3})  1 Invariants of a tensor are combinations of tensor components that do not change with changing coordinate systems. The second invariant of the strain tensor is (_{ 11} _{ 22}) + (_{ 22} _{ 33}) + (_{ 11} _{ 33}) – _{ 12}^{2} – _{ 23}^{2} – _{ 13}^{2} = (e_{1} e_{2}) + (e_{2} e_{3}) + (e_{1} e_{3}) and the third invariant is determinant[_{ ij}] = e_{1} e_{2} e_{3} 2dimensional analysis It would not seem to make any sense to have a discussion of 2D analysis after discussions of 1D and 3D analysis. Nonetheless, there are two reasons for this disorderly order. First, the 2D analysis is a special case of 3D analysis, so the discussion of 2D analysis logically follows from prior discussions. Second, we are going to use what we learn about 2D analysis to solve for the horizontal infinitesimal strain rate between three GPS sites. It makes sense to focus a bit of extra time and energy to understand the elements of the 2D analysis. Several firstorder assumptions are made in this analysis. • We assume planestrain deformation in the horizontal plane. • We assume that neglecting the vertical velocities does not significantly affect the physical interpretation of our answer, recognizing that this might be a particularly poor assumption in areas characterized by significant uplift or subsidence. • We assume that strain is homogeneously distributed across the triangular area between the three GPS sites. Our input data are locations of GPS sites and their velocities, presented in units of length per time where “time” is usually one year. The input velocity data are (effectively) instantaneous velocities, so the output data should properly be expressed as rates. In the limit as time shrinks to zero, the velocity vector and the displacement vector converge and become identical as defined in the same reference frame. Hence, we retain the usual terminology of a displacement gradient tensor rather than a velocity gradient tensor in the discussion that follows. Taking one step back, the GPS data acquisition and analysis process resolves very small displacements. The processed data available from the PBO site are daily averages of displacements relative to the initial position of the GPS receiver as resolved in a particular reference frame, such as the Stable North American Reference Frame (SNARF), relative to a particular datum, such as WGS 84. The mean velocities are computed from a time series comprised of many of these daily displacements. So we could choose to use the average displacement between the same initial time and final time for all three GPS sites (in which case our input data would be actual displacements) or the average displacement divided by the time interval associated with that displacement (in which case our input data are velocities). In the first case, our output would be a measure of the deformation during that specific time interval, and in the second case our output is the mean deformation rate. Apart from that distinction, the analysis described here is the same. The 2D analysis using data from 3 GPS sites is a perfectly constrained problem, with just the right number of known quantities to solve for the unknown quantities exactly. Of course, all of the input data have uncertainties associated with them, as does our “exact” solution. As a firstlook at infinitesimal crustal strain, the 3site trianglestrain problem provides a good foundation to build on as we learn how to work with more data. The perfectly constrained trianglestrain problem by itself is not sufficient to provide a detailed characterization of tectonic strains in an area. A full understanding of presentday crustal strain in a given area requires information from all locally available GPS sites as well as other geological/geophysical data. Due consideration would have to be paid to how a variety of nontectonic factors might affect the GPS velocity data, including flux of groundwater and surface water, igneous activity, deformation related to the earthquake cycle, mass movement that might affect a given GPS site, and so on. After we acquire the needed information from each of the GPS sites, we will know the horizontal coordinates of the initial site location (x_{o}, y_{o}) as well as the eastwest instantaneous velocity (v_{x}) and the northsouth instantaneous velocity (v_{y}) of each site (Figure 2). The site location coordinates expressed in the Universal Transverse Mercator (UTM) system are in meter units, so the instantaneous velocities are expressed in meter/year units. We will have to compute the unknown quantities: the coordinates of the translation velocity vector (t_{x}, t_{y}), the magnitude of the angular velocity vector () in nanoradian per year, and the elements of the strain rate tensor (_{xx}, _{xy}, _{yx}, _{yy}). Strain is a dimensionless quantity, but it has seemed odd to express strain rate as some number "per year." That seems to beg the question, "What, per year." So to keep us all happy and well adjusted, we throw in the word "strain" per year, even though "strain" has no dimension. The units used in infinitesimal strainrate studies using GPS data are nanostrains per year. Figure 2. Map of three GPS sites separated from one another by tens of kilometers. The eastwest and northsouth velocities of each site relative to the Stable North American Reference Frame are indicated by arrows. The velocities measured at the three GPS sites are the result of three component velocities: translation, rotation and distortion of the crust. deformation rate = translation rate + rotation rate + distortion rate The translation velocity vector has the same magnitude and direction at all three sites as the triangle of stations simply move in space relative to the reference frame. You can think of the translation velocity vector as extending from the centroid of the triangle of initial GPS site locations to the centroid of the deformed triangle. (Remember that the centroid of a triangle is the intersection of three line segments, each of which connects one apex with the midpoint of the opposite side.) The rotational velocity vector (also known as an axial vector or angular velocity vector) in this 2D case is a vertical vector extending from the center of the triangle. A positive rotation is a counterclockwise rotation. The meaning of this rotation would be clear enough if we were thinking about a rigid body that is not distorting as it is rotating. In our case, however, virtually all of the material lines in our deforming body (including the sides of the triangle between GPS sites) will be affected by shear strain, which changes the trend of the line as the area is changing shape. How do we separate the rotation of the whole area from rotations resulting from shear strain? Imagine that you use the initial locations of our GPS sites, find the centroid of the triangle, and draw a (really) big circle around that center (Fig. 3A). Then the GPS sites move, and the circle is distorted into an ellipse (Fig. 3B,D). You draw a blue line along the longer (major) axis of the ellipse, and a red line along the shorter (minor) axis that is 90° from the trend of the major axis. You use your field compass to measure the trend of the major axis in the deformed state. Then you hit the REVERSE button on any convenient time machine (remember, I started this with the word “imagine”), and reverse the deformation. You would find that the two lines you drew along the axes of the ellipse would still be 90° apart in the undeformed circle (Fig. 3C). You Figure 3. Defining rotation. A. Initial configuration of GPS triangle, with circle around the centroid of the triangle. B. Deformed triangle with circle changed to an ellipse. Major axis of the ellipse trends 60° and is blue; minor axis is red. C. Retrodeformed triangle showing orientation of red and blue material lines prior to deformation. Trend of blue line is 60°, so deformation shown in B did not involve rotation. D. Second example of deformed triangle. Major axis trends 45° and the same material line in the original triangle trended 60°, so this deformed triangle has undergone a positive (counterclockwise) rotation. might find that the blue line in the circle has the same trend as it had as the major axis of the ellipse, in which case there was no bulk rotation of the triangle during deformation (Fig. 3B). If the trend of the blue line changed during deformation, that change in angle is the bulk rotation (Fig. 3D). Distortion includes a change in shape (strain) as well as any change in volume/area (dilation). Like its counterpart described for 3D deformation, the 2D displacement gradient tensor is an asymmetric matrix that can be decomposed into the sum of a symmetric tensor and an antisymmetric tensor, representing rates of both distortion and rotation. The rotationalvelocity matrix is antisymmetric. The infinitesimal strainrate tensor _{ij} is symmetric, so that _{xy }= _{yx} and the 2D strainrate tensor thus contains only three independent variables: _{xx}, _{xy }= _{yx}, and _{yy}. Adding these two tensors together gives us the displacementrate gradient tensor e_{ij} to which we can add the translation terms {t_{x}, t_{y}} so that we account for all elements of the deformation . Note that the translation is independent of position. The complete matrix equation for the deformation at an individual site, given its initial location {x_{o}, y_{o}} and the corresponding east velocity (v_{x}) and north velocity (v_{y}) is which can be unpacked to yield the following two equations with six unknowns: t_{x}, t_{y}, , _{xx}, _{xy} and _{yy }: v_{x} = (x_{o} _{xx}) + (y_{o} _{xy}) – (y_{o} ) + (t_{x}) v_{y} = (x_{o} _{xy}) + (x_{o} ) + (y_{o} _{yy}) + (t_{y}) GPS data give us the horizontal velocities (v_{x} and v_{y}) for all three sites. Thus, we have six equations (2 for each site) with six unknowns, which makes this a perfectly constrained problem that will yield an exact answer. The six equations can then be rearranged into a single matrix equation of the form d = G m where d (for data) is the matrix of the known velocity components for our three stations, m (for model) is the matrix of our six unknowns, and G (named after George Green) is the matrix of coefficients that relates d to m using the equations that we just developed. Because we already know all of the elements of the data matrix d, but do not know the model parameters in matrix m, we will need to solve the inverse problem so that all of the known quantities are collected on one side and the unknowns are on the other side of the equation. m = G^{1}d Once we have computed matrix m, its six components include the coordinates of the translational velocity vector (in m/yr), the rotational velocity (in radian/yr), and the three independent elements of the strain matrix. The principal strain axes are found by computing the eigenvectors of the 2D strain tensor _{ij}. The eigenvectors are unit vectors in the directions of the principal strain axes. The eigenvalues of _{ij} are the principal extensions in the principal directions. The larger of the two eigenvalues is the greater principal extension, e_{1}, and the length of the semimajor axis of the horizontal strain ellipse is equal to the greater principal strain, S_{1} = e_{1} + 1. The semiminor axis is S_{2} = e_{2} + 1, where e_{2} is the lesser principal extension. In the specific case in which we have a 2D strain tensor the eigenvalues are given by . The larger eigenvalue is _{1} (also known as e_{1}) and the smaller eigenvalue is _{2} (also known as e_{2}). An eigenvector corresponding to _{1} is , and an eigenvector associated with _{2} is . The unit eigenvectors can then be determined by dividing each of the components of these vectors by their length or norm. For example, the unit eigenvector associated with eigenvalue _{1} is The magnitude of the infinitesimal shear strain (_{max}) at 45° to the maximum principal strain axis can be computed a couple of different ways, given the data we now have in hand. _{max} = _{max} = e_{1}  e_{2} The area strain, which is also the first invariant of the 2D strain tensor, is given by area strain = _{11} + _{22} = e_{1} + e_{2} = (S_{1} S_{2})  1 The second invariant of the 2D strain tensor is (_{ 11} _{ 22}) – _{ 12}^{2} = e_{1} e_{2} and the third invariant is determinant[_{ ij}] = e_{1} e_{2} Remember that the invariants of a tensor are combinations of tensor components that do not change with changing coordinate systems. Estimating the uncertainty of our strain solution The uncertainty of our strain estimate is a function of the individual uncertainties in the velocity estimates (which are available from the same source as the GPS data) and how these uncertainties map into the model space. The data uncertainty can be expressed as a covariance matrix (covd) with diagonal values equal to the individual variances (^{2}), assuming that the error in each station’s velocity estimate is independent of the other stations. For a linear equation of the form m = G^{1}d (our inverse solution). The covariance of the model (cov m) can then be expressed as (Menke, 1989) or more simply as where the inverse of the data covariance matrix (cov d)^{1} is simply a diagonal matrix with values equal to 1/^{2} for each of the velocity estimates. The diagonal elements of the 6x6 model covariance matrix (cov m) are the variances (^{2}) for each of the model parameters we have estimated. Their standard deviation () can be approximated by taking the square root of variance. (Note that the offdiagonal elements of the model covariance matrix (cov m) might not be zero, indicating some covariance between the various model parameters.) The diagonal values in rows 1 and 2 of matrix (cov m) are the variances for the x and y coordinates, respectively, of the translation vector. The diagonal value in row 3 is the variance of the angular velocity vector, . The diagonal values in rows 46 are the variances for _{xx}, _{xy} and _{yy}, respectively. Now we are ready to calculate the translation, rotation, and strain rate for a set of GPS stations as well as estimate the uncertainty of our results. References Allmendinger, R.W., Cardozo, N, and Fisher, D.M., 2012, Structural geology algorithms, vectors and tensors: Cambridge University Press, 289 p., ISBN 9781107401389. Cardozo, N., and Allmendinger, R.W., 2009, SSPX  A program to compute strain from displacement/velocity data: Computers and Geosciences, v. 35, p. 13431357, doi: 10.1016/j.cageo.2008.05.008. Jaeger, J.C., 1964, Elasticity, fracture and flow: Methuen, New York, 212 p. Malvern, L.E., 1969, Introduction to the mechanics of a continuous medium: Englewood Cliffs, New Jersey, PrenticeHall, Inc., 713 p. Means. W.D., 1976, Stress and strain – Basic concepts of continuum mechanics for geologists: Englewood Cliffs, New Jersey, PrenticeHall, Inc., 339 p. Menke, W., 1984, Geophysical data analysis – Discrete inverse theory: Orlando, Florida, Academic Press, 260 p. Menke, W., 2012, Geophysical data analysis – Discrete inverse theory MATLAB edition: Orlando, Florida, Academic Press, 348 p. Oertel, G., 1996, Stress and deformation – A handbook on tensors in geology: New York, Oxford University Press, 292 p., ISBN 0195095030. Pollard, D.D., and Fletcher, R.C., 2005, Fundamentals of structural geology: Cambridge University Press, 500 p. Savage, J.C., Gan, W., and Svarc, J.L., 2001, Strain accumulation and rotation in the Eastern California Shear Zone: Journal of Geophysical Research, v. 106, B10, P. 21,99522,007. Timoshenko, S.P., and Goodier, J.N., 1970, Theory of elasticity [3^{rd} edition]: New York, McGrawHill Book Company, 567 p., ISBN 070647208. Version of 6 September 2012 Page 