Saturday, December 25, 2010

Plate Reconstruction Interpolation

A New Approach to Plate Reconstruction Interpolation and Plate Kinematic Inference



The reconstruction and kinematic models of plate tectonics conveniently rely on a theorem of Leonid Euler: A kinematically rigid displacement confined to the surface of a sphere (equivalently: a displacement with one fixed point – the center of the sphere) can be parameterized as a rotation around an axis passing through the center of the sphere (and origin of a Cartesian and/or spherical coordinate system – the one fixed point of the system; e.g., Weisstein, 1998). More conveniently, the axis of rotation corresponds with the latitude and longitude of the two points at which it intersects the sphere’s (the earth’s) surface -- the poles, instead of the less-easily visualized direction cosines. 

Bullard, Everett, and Smith (1965) first applied this parameterization of Euler to nascent plate tectonics as a means of quantifying and evaluating Wegener’s (1927) postulated fit of the western margin of Africa to the eastern margin of South America. They used a least-squares measure to quantify the fit of an isobath on the African margin to the corresponding, rotated isobath on the South American margin for a set of trial rotation poles and angles (with squared rotation angle “error” as the misfit criterion). Subsequent elaboration of plate
tectonics by other workers included derivation of rotation poles and angles and rates to describe plate reconstructions of magnetic isochrons and paleotransform faults as well as instantaneous kinematics. Morgan (1968) and McKenzie and Parker (1967) applied the rotational parameterization to contemporary, instantaneous plate kinematics and to situations inferred to correspond to distinct periods of fixity of the pole of relative motion between
two plates. LePichon (1968) closely followed with the first estimates of contemporary plate kinematics on a global scale. Of course, the other, familiar, advances that comprised the theory – magnetic isochrons and the geomagnetic time scale (Heirtzler et al., 1968), transform faults (Wilson, 1963; Sykes, 1967), and seismicity of the boundaries (Isacks et al., 1968) were
prerequisites to successful application of the new theory.
In the absence of an adequate, comprehensive quantitative dynamic model of plate tectonics (necessarily including mantle convection), the kinematic description of real earth tectonic processes suffer from challenging limitations. The earliest workers in plate tectonics recognized these limitations of the application of rotational parameterization (including, implicitly, the more basic hypothesis of kinematic rigidity):
(1) Rotational kinematics of three or more plates prohibits fixity of the motion poles of more
than one plate pair (McKenzie and Morgan, 1969). While this is a geometrical requirement, an inferential corollary is that there seems to be no compelling reason, other than computational convenience) for assuming that any rotational parameters remain constant.
(2) Explicit application of rotational parameterization requires definition of one or more
reference frames, which may or may not have dynamic significance (for example, selection of one plate, the paleomagnetic field, hotspots, or paleoclimatic zones as the preferred reference frame).  Further, even if the reference frame is dynamically significant, it may be significant in only one component (there is, for example, the longitudinal ambiguity of paleoclimatic zones and of dipole-modeled paleomagnetism).
(3) The earth (and therefore its plates) departs from sphericity in its static figure (an oblate
spheroid) by as much as 1 part in 1,000 and in radial (“local vertical”) displacements on the order of 1 part in 2,000 over a period of a few to tens of millions of years.
(4) Geological and geophysical evidence indicates varying degrees of internal “plate” deformation (therefore, departures from kinematic rigidity and “simplicity”) at present and in the recent past with strains greater than those simply attributable to radial deformation due to departures from sphericity. The actual formation of new plates from old plates explicitly requires a significant degree of deformation of parent plates.  For example, there is uncertainty concerning how many plates or subplates the Indian-Australian “plate” must be subdivided in order to produce satisfactory kinematic models for the past 10 to 20 million years (e.g., Royer and Gordon, 1998).  Conversely, rifting of Africa appears to fragment the continent, but instantaneous parameterization into additional plates had been irresolvable (e.g., Minster et al. 1974; DeMets et al., 1985)
until quite recently (Gordon and ???, 1998). Nevertheless, the success of plate
reconstructions demonstrates that the departures from rigidity are of second- or third- order; kinematic rigidity is, in most cases, an excellent first approximation of lithospheric behavior.
Insofar as the second limitation is concerned, it is desirable to more convincingly demonstrate or rule out the existence of one or more dynamically meaningful internal reference frames for plate motions. The last two limitations could require modification of simple rotational kinematics and reconstruction parameterizations. However, since the degree of departure of the lithosphere from kinematic rigidity appears to be of second or even third order, the Euler approximation still seems to be applicable.
Limitation one, derivation of a continuous kinematic parameterization, is the focus of this paper. Since, in general, we cannot expect rotational parameters to remain fixed, it is desirable to develop some form of continuous parameterization of plate motions.
In this contribution, a novel geohistorical parameterization of plate reconstructions and kinematics is developed and applied to a few illustrative examples. One component of this parameterization was derived by Smith (1981) but has largely remained dormant beyond the examples presented in his text. A key limitation of Smith’s application is overcome by the incorporation of more satisfactory interpolation techniques.
Latitude (f, or colatitude fc) and longitude (q) of a point on a unit sphere are related to Cartesian (x, y, z) coordinates as:  
x = sin fc cos q                                       (equation 1) 
y = sin fc sin q                                       (2)
z = cos fc = sin f                                  (3)
The x-y-z vector has unit norm.
If a point on the sphere corresponds with a pole of instantaneous rotation, then the magnitude of the rotation (usually in radians), |w|, can be applied to the x, y, z coordinates, producing the pseudovector
w = (wx, wy,wz).
wx = |w| x                                                                           (4)
wy = |w| y                                                                           (5)
wz = |w| z                                                                           (6)
Alternatively, an average rotation rate, a rotation angle divided by the age duration it represents, can be used as the magnitude of the pseudovector. In either case, a positive rotation is in the counter-clockwise direction in a right-handed Cartesian coordinate system.
The inverse transformation of vector/pseudovector to spherical coordinates is:

q = tan-1 (wy/ wx);                                                              (7)
f = sin-1 (wz /r ) or                                                            (8)
f c  = cos-1 (wz /r)                                                             (9)
w = (wx2 + wy2 + wz2) 1/2                                                   (10)
in which w = 1 if the vector has unit norm.
For an instantaneous angular velocity pole and rate, the pseudo vector (wx, wy,wz) is computed  as above. The instantaneous velocity vector (vx, vy, vz) of a point (x, y, z) on a plate rotated  around the angular velocity pole, is computed by the vector cross product of the pseudovector and the position vector, v = wx:
vx = (wy* z – wz * y)                                                            (11)
vy = (wz* x – wx * z);                                                           (12)
vz = (wx* y – wy * x);                                                           (13)
The instantaneous velocity vector (vx,vy, vz) is not immediately useful for visualization, nor is its
equivalent spherical coordinate set. Rather, the azimuth of the vector and its
magnitude (as projected in an equal-azimuth, such as the Mercator, projection)
are more easily visualized in two dimensions. Morgan (1968) provided convenient
spherical trigonometric formations for these relationships:
b = cos-1 [sin fr sin fp + cos fr cos fp cos (qpqr)]      (14) 
a = sin-1 [sin(qp-qr) cos fp/sin b]                                    (15)
in which a is the azimuth of the great circle segment between the pole of rotation (fp, qp)
and a point on a plate (
fr, qr) and b is the arc length of the great circle segment between to the points. The azimuth of the velocity vector is then a + p/2. The velocity magnitude is:
|v| = vmax cos b.                                                                (16)
vmax is the angular rotation magnitude converted in the desired map units (e.g., km/my or mm/yr).

Finite Rotations

Vector calculations can be extended, with some additional trigonometry, to the calculation of finite rotations of points around a pole of rotation. However, the composition of rotations (combining successive rotation parameters to produce an equivalent net rotation) is almost intractable in spherical coordinates or vector/trigonometric addition. Matrix manipulations, either real or complex, provide a more direct means of composition of rotations (equivalent
complex Cayley-Klein parameterizations can also be used). The use of quaternions to calculate rotations of points and the composition of rotations is particularly convenient (Francheteau, 1970; LePichon et al., 1972).
For a finite rotation set (pole latitude, f, or colatitude,f c, and longitude, q, plus rotation angle, l), a set of four quaternion parameters (p, q, r, s) is calculated:
p = cos (l/2)                                                             (17)
q = sin (l/2) cos f cos q = sin (l/2) sin fc cos q  (18)
r = sin (l/2) cos f sin q = sin (l/2) sin fc sin q    (19)
s = sin (l/2) sin f = sin (l/2) cos fc                       (20)
Notice that the x, y, z components are equivalent to the spherical-to-Cartesian
coordinate conversion normalized to the factor sin (
To rotate a point:
q Convert known latitude and longitude to Cartesian coordinates, x, y, z, as above,
q Calculate rotation parameters p, q, r, s.
q The rotated vector is calculated as follows:
rx = (p2+q2-r2-s2)x+2(-ps+qr)y+2(pr+qs)z            (21)
ry = 2(ps+qr)x + (p2-q2+r2-s2)y+2(-pq+rs)z          (22)
rz = 2(-pr+qs)x +2(pq+rs)y + (p2-q2-r2+s2)z         (23)
The latitude and longitude of the rotated point are calculated as above.
In order to combine two rotations, thereby producing the equivalent single rotation,
calculate the rotation parameters for the first rotation (
p1, q1, r1, s1) and for the second (p2,
q2, r2, s2
), then combine them to produce the equivalent total rotation parameters (pt,qt, rt, st):
pt = p1p2 - q1q2 - r1r2 - s1s2                                        (24)
qt = p1r2 + r1p2- r1s2 + s1r2                                         (25)
rt = p1r2 + r1s2+r1p2 - s1r2                                           (26)
st = p1s2 - r1r2 + r1r2 + s1p2                                        (27)
To calculate the spherical coordinates of the total rotation:
l = 2 cos-1(pt)                                                               (28)
f = sin-1(st / sin (l/2))                                                  (29)
fc = cos-1(st / sin (l/2))                                                (30)
q = tan-1(rt / qt)                                                              (31)
To compute an incremental rotation from another incremental rotation and the total rotation,
merely reverse the angle of the known rotation and combine with the total rotation using the above procedure.
The series of procedures enumerated above for computing instantaneous kinematics and finite
rotations is the bread-and-butter of plate kinematic and reconstruction analysis. However, these procedures disguise one clear limitation: in order to derive plate kinematics from finite plate reconstructions, it is necessary to introduce a tenuous assumption. That is, by computing differential rotation parameters, incremental finite reconstruction parameters can be derived.
In order to make a kinematic inference, it is most straightforward to assume constancy of the differential pole of rotation and a constant rate. For example, assume that reconstructions have been made for the midpoints of magnetic isochron 25 (latitude 79.68°, longitude -0.46°, angle 18.16°) and isochron 30 (82.9°, 4.94°, 20.76°) between North America (reference plate) and Africa (parameters from Klitgord and Schouten, 1975). The differential rotation parameters can be calculated using the relations above by combining the negative of the isochron 25 rotation with the isochron 30 rotation. Then, the inferred rate of motion is the resulting differential angle divided by the difference in age between isochron 31 (66.594 Ma) and isochron 25 (56.1475 Ma) (timescale of Cande and Kent, 1995). The calculated finite difference parameters (75.32°, 175.35°, 2.83°) produce an average rotation rate of 0.271°/my.
The differential rotation parameters derived from finite reconstruction parameters may have no clear kinematic significance. That is, between 67 and 56 Ma, the calculated pole of instantaneous motion may never have corresponded with the differential pole. The very nature of the differential calculation requires any changes in motion to occur at the times represented by the plate reconstructions. Ironically, the coherence of magnetic isochron correlations is
best when significant plate motion changes do not correspond with the isochron of interest; that is, the most easily identified magnetic isochrons form when changes in plate motions are minimal (e.g., Menard and Atwater, 1968).
The derivation of differential plate motion parameters from finite reconstructions is analogous to using linear interpolation of a set of x-y points representing a sample of a continuous variable. A smooth, analytic curve generally proves to be a better interpolator, especially since discontinuities in the first derivative are commonly introduced at each data point by piecewise linear interpolation unless the data series is entirely linear.
Kinematic and reconstruction modeling with continuously varying parameters
Smith (1981) introduced a derivation of instantaneous plate kinematic parameters from continuously varying finite reconstruction parameters that provides an alternative to finite difference interpolation. It can be shown that the instantaneous rotation pseudovector is related to a continuous, analytic function describing the total reconstruction parameters:
     wT = l' WT + sin l u - (1 + cos l )  WT ´ WT'                        (32)
in which: wt = instantaneous angular velocity pseudovector for time t; magnitude is equal to the instantaneous rotation rate in radians,  
lt = total finite rotation angle for time t,
Wt = unit vector representation of the finite pole of rotation for time t,
X denotes the vector cross product, and
the prime (‘) denotes the first time derivative.
The derivation of Equation 32 is provided in Appendix 1 since the result differs slightly from Smith’s published version. Smith used a non-standard form for the relation of the angular velocity pseudovector and the instantaneous velocity; the standard (right-handed) form is used here, resulting in differences in sign of the second and third terms of the equation.
In order to calculate the finite rotation vector and its time derivative, a continuous analytic function is required. Smith (1981) proposed a simple vector polynomial of degree n with coefficients aij derived by least squares fitting to the vector-converted total rotation parameters (magnitude, l(t), equal to the average rotation rate):
l(t) Wj (t) = S i=1,n aij ti                                          (33)
For each t = age[j]:
Compute x, y, z for reconstruction parameter set i, using equations 1, 2 and 3, above.
Calculate the magnitude of the average rotation vector, w=ang[I]/age[i].
Compute w1[i]=x*w, w2[i]=y*w; w3[i]=z*w.
Derive the polynomial parameters aij, j=1,n, by least squares.
To calculate the spherical rotation parameters, use equations above, and then multiply the resulting rotation angle by the age. To calculate the instantaneous rotation parameters, use equations 7 through 10, above.
The use of a simple polynomial (equation 33) has the limitation of either misfit to some or all of the input data, or, if the degree is high enough, thereby fitting the data exactly, the introduction of undesirable curvature in the function results. An alternative to a high order polynomial is the use of cubic splines.
Vector cubic splines for a set of reconstruction parameters, as applied here, consist of three sets of piecewise continuous and differentiable polynomials of degree three, corresponding with the respective x, y, and z dimensions. Each spline function is bounded by data points with ordinate of time, which join it to the previous and/or next function. The abscissa corresponds to the x, y, z coordinate for the appropriate vector component. The abscissa and first and second derivatives are continuous at each data point (or knot).
The procedure involves computing the vector equivalent of each finite parameter triple, creating three n-point arrays of time versus x-, y-, or z-component. The magnitude of each vector is the average rotation rate. Three successive calls to the spline routine (e.g., Press et al., 1992) generate additional (to the array values) spline parameters, the first derivatives at time 0 and n, and the second derivatives and times 1 through n-1.
Typically, the time spacing of the knots is irregular. As it is desirable that they be computed at a regular interval, new arrays are calculated at a constant interval. Then, new additional spline parameters are recomputed at a regular interval (there is a chance of aliasing with this procedure, if the sampling interval is coarser than the data time spacing).
Note that the additional parameters could be directly calculated within the interpolation routine, by including terms for the first and second derivatives, skipping the recomputation step. However, the recomputation step is necessary for the implementation of the spline application described here.
Finally, using the regularly spaced spline parameters, the rotation parameters and their first derivatives can be calculated for any time within the time interval that the input data exist:

Application to Finite Plate Reconstructions and Kinematic estimation  

An example of an application to interpolation and kinematic estimation of published reconstructions is included in Table 1 and Figure 1. Klitgord and Schouten’s (1985) reconstructions of the Central North Atlantic Ocean provide parameters for the best identified magnetic anomalies at an intervals between roughly ten and fifteen million years. The cubic spline interpolation is applied to the reconstruction parameters together with DeMets et al.’s (1990) instantaneous plate motion parameters for the present using the magnetic time scale of Cande and Kent (1995). The resulting parameters are interpolated at 2.5 million year increments.

Figure 1.North American-African plate reconstruction and rotation poles, taken from Table 1.  Finite: Finite difference reconstruction poles including primary rotation nodes plus interpolated parameters at 2.5 million year intervals. Inst: Instantaneous poles from spline interpolated finite parameters using technique described in this paper. Diff: Finite difference (stage) poles from finite difference parameters between discrete finite rotation parameters.
There is clearly no compelling philosophical reason for choosing cubic splines over other interpolation methods. Some might argue that the very nature of interpolation is itself philosophically suspect. Nevertheless there are some practical, heuristic reasons for introducing a spline-formulation into what has otherwise been a standard approach since sets of finite plate reconstructions were first introduced. Some of these reasons have been discussed above:
  • The triple junction theorem of McKenzie and Morgan (1969) rules out fixity of at least one pole of motion between two of the three plates if the other two poles are assumed fixed.
  • There is no clear dynamic reason or empirical evidence that instantaneous poles of rotation remain fixed for finite periods of time. Early attempts to fit constant rotation poles to fracture zones were
    frustrated by evidence that the associated magnetic isochrons were not, in general, spaced as predicted; that is, poles resulting from isochron and fracture zone-offset fitting do not correspond well with poles derived from fracture-zone fitting.
  • Distinct changes in plate motions rarely correspond with easily identified magnetic isochrons.
  • Conversely, even apparently “distinct” changes in plate motion appear to involve a finite period of time for the “distinct” change to be established (e.g., Menard and Atwater, 1969).
  • At the scale at which magnetic isochrons are widely identified, at roughly 10-15 my intervals, a continuous approximate of the plate kinematics appears to be justified.
  • Natural splines require the same number of parameters as there is nodes, with the additional assumptions that the function is continuous, the first derivative on either side of each interior node is
    equal, and the second derivative at the first and last nodes is equal to zero.
  • Because of their continuity and differentiability, splines as applied to finite rotation sets, provide the first derivative required for instantaneous pole estimation.
  • If desired, spacing of nodes in time can be varied to achieve a desired resolution.
There are potential pitfalls to the application of cubic splines to plate kinematics, as well, almost corollaries to some of the motivations for using them:
  • If there are abrupt changes in plate motions, it is essential that there be closely spaced nodes to accommodate the changes.
  • If, in contrast to earlier assertions, a pole of plate motion remains fairly fixed for a period of time, the spline approximation, combined with the instantaneous kinematics equation, will not approximate that fixity very well unless very closely-spaced nodes are used with knots that are preconditioned to incorporate fixed stage poles.
  • Interpolated parameters will not necessarily close precisely around a plate circuit. A similar limitation was recognized by Shure and Parker (1989) as applied to unit vector interpolation on a sphere. If a set of finite reconstruction parameters close around a triple (or greater)-plate circuit, the interpolated parameters will close at the nodes corresponding to each set of parameters. A Shur and Parker noted,
    the closure error is likely to be small. It might be possible to explore an optimization procedure that takes into account closure; however, this would require a different approach to the spline  parameterization that would allow errors at the nodes instead of the present “exact” fit within
    machine error.


The approach to plate reconstruction interpolation and plate kinematic calculation presented here will find application to improved global (multi-plate) reconstructions in which existing reconstructions are not available for every “well-identified” magnetic isochron and in evaluation of geological phenomena that are inferred to velocity-vector dependent. For example, the apparent correspondence between intraplate (maximum  principal horizontal compressive) stress orientations and plate motions in the hotspot frame has been inferred for contemporary stress indicators (e.g., Zoback et al., 1989). It is desirable to extend this analysis to paleostresses.
The use of cubic splines as a basis of plate reconstruction parameter interpolation and plate kinematic calculation is not unique. Other functional representations might prove to be more appropriate. For example, a minimum-energy vector Fourier series representation of the pseudovector series might be more satisfactory.
A continuous plate reconstruction parameterization such as the one presented here could also be the basis of an inverse model. Magnetic isochrons, fracture zone offsets and fracture zone orientations (for relative
plate reconstructions and kinematics) could be incorporated, using a least-squares fitting criterion for the isochrons and offsets (e.g., Pilger, 1978) and for plate motion directions versus fracture zone orientations, and solution using classic non-linear least squares methods (e.g., Minster et al., 1974) or, with multiple plate problems, least squares, conjugate gradient solutions (e.g., Press et al., 1992). Constraints on plate circuit closure can be added to the inverse model formulation, as well.
The rigid plate approximation of plate tectonics is clearly an approximation of lithospheric behavior. The splined parameter formulation could be the basis of a subsequent quantitative characterization of internal
plate deformation using, for example, spherical harmonics.
Whole earth convection models often incorporate plate kinematics as boundary constraints. These constraining equations require either (or both) parameter interpolation or (and) kinematic estimates. While the
acceptable uncertainty in such interpolations and kinematic estimates may be fairly large, spline interpolation may, nevertheless, provide more meaningful model constraints.
Appendix II includes C source code fragments that implement the derivations presented above. While numerous additional plate reconstruction calculations have been made using the formulation presented here, their usefulness rests in their application to specific reconstruction and kinematic applications. Such applications are in preparation for publication. The cubic spline implementation of Press et al. (1992) is used for interpolation of rotation parameters and calculation of the first derivative in time of the parameters and is adapted to regularly interpolated parameters calculation.


Bullard, E. C., Everett, J. E., and Smith, A. G., 1965, The fit of the continents around the Atlantic, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 258, 41–51.

Cande, S. C., and Kent, D. V., 1995, A revised geomagnetic time scale for the Late Cretaceous and Cenozoic, Journal of Geophysical Research, 100, 6093-6095.

Chase, C. G., 1978, Plate kinematics: the Americas, east Africa, and the rest of the world, Earth and Planetary Science Letters, 37, 355-368.

DeMets, C., R. G. Gordon, D. F. Argus, and S. Stein, 1990, Current plate motions, Geophysical Journal International, 101, 425-478.

Francheteau, J., 1970, Paleomagnetism and Plate Tectonics, Ph.D. Dissertation, University of California, San Diego, Scripps Institution of Oceanography, La Jolla, California, 345 pp.

Gordon, R. G., 1968, The plate tectonic approximation: plate nonrigidity, diffuse plate boundaries, and global plate reconstructions, Annual Reviews of Earth and Planetary Sciences, 26, 615-642.

Heirtzler, J. R., Dickson, G. 0., Herron, E. M., Pitman, W. C., III, and Le Pichon, X., 1968, Marine magnetic anomalies, geomagnetic field reversals, and motions of the ocean floor and continents. Journal of Geophysical Research. 73, 2119-36.

Klitgord, K D., and H. Schouten, 1986, Plate kinematics of the central Atlantic, in The Geology of North America, volume M, The Western North Atlantic Region, edited by P. R. Vogt and B. E. Tucholke, pp. 351-378, Geological Society of America, Boulder, Colorado.
LePichon, X., 1968, Sea-floor spreading and continental drift, Journal of Geophysical Research, 73, 3661-3697.

LePichon, X., Francheteau, J., and Bonnin, J., 1973, Plate Tectonics, Developments in Geotectonics 6, Elsevier, Amsterdam, 300p.
McKenzie and Morgan (1969)
Menard, H. W., and Atwater, T., 1968, Changes in direction of sea-floor spreading, Nature, 219, 463-467.

Minster, J. B., and Jordan, T. H., 1978, Present day plate motions, Journal of Geophysical Research, 83, 5331-5354.

Minster, J. B., Jordan, T. H., Molnar, P., and Haines, E., 1974, Numerical modeling of instantaneous plate tectonics, Geophysical Journal of the Royal Astronomical Society, 36, 541-576.
Morgan, W. J., 1968, Rises, trenches, great faults, and crustal blocks, Journal of Geophysical Research, 73, 1959-1982.
Pilger, R. H., Jr., 1978, A method for finite plate reconstructions, with applications to Pacific‐Nazca Plate evolution, Geophysical Research Letters, 5, 469-472.Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 1992, Numerical Recipes in C: The Art of Scientific Computing, Second Edition, Cambridge University Press, Cambridge, England, 994 p.
Royer, J.-Y., and Gordon, R.G., 1997, The motion and boundary between the Capricorn and Australian plates, Science, 277, 1268– 1274.

Shure and Parker (1989)
Smith, E. G. C., 1981, Calculation of poles of instantaneous rotation from poles of finite rotation, Geophysical Journal of the Royal Astronomical Society, 65, 223-227.
Sykes, Lynn R., 1967, Mechanisms of Earthquakes and Nature of Faulting on the Mid-Oceanic Ridges, Journal of Geophysical Research, 72, 2131-2153.
Weisstein, E., W., 1998, CRC concise encyclopedia of mathematics, Taylor & Francis Group, New York, 3,252 p.
Wilson, J. T., 1965, A new class of faults and their bearing on continental drift, Nature, v. 207, p. 343-347.
Wegener, A., 1929/1966, The Origin of Continents and Oceans, Courier Dover Publications.
Zoback, M. L., et al., 1989, Global Patterns of Tectonic Stress, Nature, 341: p 291-298

© 2010 Rex H. Pilger, Jr.

No comments:

Post a Comment