Origin and Application of the Twinned Calcite Strain Gauge

This paper is a personal account of the origin and development of the twinned-calcite strain gauge, its experimental verification, and its relationship to stress analysis. The method allows the calculation of the three-dimensional deviatoric strain tensor based on five or more twin sets. A minimum of about 25 twin sets should provide a reasonably accurate result for the magnitude and orientation of the strain tensor. The opposite-signed strain axis orientation is the most accurately located. Where one strain axis is appreciably different from the other two, that axis is generally within about 10° of the correct value. Experiments confirm a magnitude accuracy of 1% strain over the range of 1–12% axial shortening and that samples with more than 40% negative expected values imply multiple or rotational deformations. If two deformations are at a high angle to one another, the strain calculated from the positive and negative expected values separately provides a good estimate of both deformations. Most stress analysis techniques do not provide useful magnitudes, although most provide a good estimate of the principal strain axis directions. Stress analysis based on the number of twin sets per grain provides a better than order-of-magnitude approximation to the differential stress magnitude in a constant strain rate experiment.


Introduction
This paper is a personal account of the origin and development of the calcite strain gauge calculation technique along with a discussion of subsequent experimental studies that tested and added to the concept. Key aspects of the original work will be evaluated and the method, compared to contemporaneous stress analysis techniques.
I developed the twinned-calcite strain gauge calculation technique while studying for my PhD with Bill Chapple at Brown University. Bill was one of the few academic structural geologists then applying mechanics to interpret structures and his classes introduced us to stress and strain models and rock mechanics, as well as to finite-amplitude buckle folding theory. I decided to measure strain in natural buckle folds to see how it fit with the mathematical models. This was to be conducted using a measurement technique probably known nowhere else, the method of Conel [1]. James Conel completed his PhD at Cal Tech under the supervision of Barkley Kamb, and Bill Chapple was also a Kamb student. Conel developed a form of a calcite twin strain gauge in his dissertation but never published it. Fortunately, Bill had a copy of the dissertation because there was almost no other way that I or anyone else was likely to have learned about it. A bit of serendipity is that John Spang was already at Brown studying for his PhD with Chapple. John had used a U-stage in his master's research; thus, he knew how to operate the instrument, and he taught me how to do it. Together, we learned how to measure calcite twins (later the two of us were occasionally known as the calcite twins). Another bit of serendipity was that John Imbrie was in the department teaching a class on advanced statistical methods.
This paper approximately follows the chronology of the development of the technique. It begins with the conceptual background that gave rise to the strain gauge technique, then its development and the results of the initial testing using synthetic and field data. This is followed by the results of the first test of the method using experimentally deformed

Discovery and First Results 1968-1972
In the late 1960s, Turner's 1953 dynamic analysis [2] of calcite based on e-twins was well known. Turner's dynamic analysis uses the crystallographic orientation of calcite twins to find the compression (C) and tension (T) axes of stress most favorably oriented to cause the twinning (Figure 1). Handin and Griggs [3] had previously shown that the amount of shear strain in a completely twinned grain could be computed from the crystallography. Conel [1] was the first to recognize that calcite twins in a partially twinned grain represent a measurable simple shear strain ( Figure 1). The amount of strain recorded by a twin set is a function of the number of twins, their widths, and the width of the grain in the direction normal to the twin plane. Each twin set is interpreted to represent an independent strain gauge embedded in the deforming rock. The engineering shear strain magnitude, the measure most commonly used by geologists, for a single twin set is tan ψ (Figure 1), and the tensor shear strain is 1 /2 the engineering shear strain [4].
data. This is followed by the results of the first test of the method using experimentally deformed limestone. Then, the issues of recognizing multiple deformations and the minimum number of measurements for a reliable result are considered. The paper concludes with discussions of stress analysis based on new compilations of the original experimental data, strain partitioning in a natural fold, and some suggestions for future research.

Discovery and First Results 1968-1972
In the late 1960s, Turner's 1953 dynamic analysis [2] of calcite based on e-twins was well known. Turner's dynamic analysis uses the crystallographic orientation of calcite twins to find the compression (C) and tension (T) axes of stress most favorably oriented to cause the twinning (Figure 1). Handin and Griggs [3] had previously shown that the amount of shear strain in a completely twinned grain could be computed from the crystallography. Conel [1] was the first to recognize that calcite twins in a partially twinned grain represent a measurable simple shear strain ( Figure 1). The amount of strain recorded by a twin set is a function of the number of twins, their widths, and the width of the grain in the direction normal to the twin plane. Each twin set is interpreted to represent an independent strain gauge embedded in the deforming rock. The engineering shear strain magnitude, the measure most commonly used by geologists, for a single twin set is tan ψ (Figure 1), and the tensor shear strain is ½ the engineering shear strain [4].
Once the strain is measured in a number of grains, the next question is how to combine this information to find the bulk strain tensor. Conel did a 3D tensor transformation on the strain in each twinned grain to bring the strain magnitudes all into the same (thin section) coordinate system. Then, he averaged each of the six tensor components, weighting the strain from each grain by the area of the grain. This produced strain tensors with magnitudes and directions that looked reasonable in the fold he studied.  [5]). e is the pole to the twin plane and +g is the glide direction: C = compression axis, T = tension axis, ti = thickness of a twin, pi = offset of boundary caused by a twin, q = total boundary offset due to twinning, α = rotation angle from untwinned to twinned position (38°17′), ψ = angle of shear, Γeg = tensor shear strain in the +e and +g coordinate system. Conel did the tensor transformation and averaging calculations by hand. When I began to use his strain technique it was then possible to use a computer. Researchers at Brown all had ready access to the university mainframe computer. To do a calculation, I wrote a program in FORTRAN, punched it into cards, submitted the cards, and then a day or so later the results, and the cards came back. Debugging the program was a very slow process. To validate the program, I decided to generate a synthetic data set. I produced a set of random twin orientations and, from an initial strain tensor, calculated the strain for each twin set in the deformed state. This was an error-free data set for which a correct  [5]). e is the pole to the twin plane and +g is the glide direction: C = compression axis, T = tension axis, t i = thickness of a twin, p i = offset of boundary caused by a twin, q = total boundary offset due to twinning, α = rotation angle from untwinned to twinned position (38 • 17 ), ψ = angle of shear, Γ eg = tensor shear strain in the +e and +g coordinate system.
Once the strain is measured in a number of grains, the next question is how to combine this information to find the bulk strain tensor. Conel did a 3D tensor transformation on the strain in each twinned grain to bring the strain magnitudes all into the same (thin section) coordinate system. Then, he averaged each of the six tensor components, weighting the strain from each grain by the area of the grain. This produced strain tensors with magnitudes and directions that looked reasonable in the fold he studied.
Conel did the tensor transformation and averaging calculations by hand. When I began to use his strain technique it was then possible to use a computer. Researchers at Brown all had ready access to the university mainframe computer. To do a calculation, I wrote a program in FORTRAN, punched it into cards, submitted the cards, and then a day or so later the results, and the cards came back. Debugging the program was a very slow process. To validate the program, I decided to generate a synthetic data set. I produced a set of random twin orientations and, from an initial strain tensor, calculated the strain for each twin set in the deformed state. This was an error-free data set for which a correct calculation program would return the original strain magnitudes and directions. However, Conel's calculation technique did not give the starting strain magnitudes or the orientations of the principal axes. Was this a programming error or something different?
The equation for the strain in a grain of a given orientation, Γ eg , in the 3D tensor transformation equation has six components [4]. After a week or two of perplexity, I realized that this represented one equation in six unknowns, which could be solved with measurements from multiple grains. As twinning alone does not cause or measure volume change, one unknown is eliminated by assuming constant volume. Thus, the equation could be solved uniquely with five equations from five differently oriented twin sets. With more twin sets, a best-fit solution can be found by multiple linear regression (least squares, for short), along with a statistical measure of the goodness of fit. The least-squares method finds the best-fitting strain ellipsoid for the twin-set strains in their measured U-stage orientations. The eigenvalues and eigenvectors of the result are the principal strain directions and magnitudes. The original strain tensor was obtained from the synthetic data by this calculation, proving that the least-squares approach worked and validating the computer program. This is the origin of the calcite strain gauge, the need to debug a complex computer program. (In 1994, Mark Evans brought the program into the desktop computer era and added the ability to plot the directional information on stereonets [6]; see Tielke et al. in this issue for a new Excel version).
Forward modeling of randomly oriented potential twin sets produced grains with negative shear strain. This is exactly opposite to the shear direction shown in Figure 1 and is impossible. By definition, the strain in a twin set can only be positive. In the initial modeling, grains with negative twin strains were simply discarded, but negative twin strains turned up later in the field data and would need to be explained.
Having error-free synthetic twinned grains made it possible to compare the different methods for interpreting the data. Comparisons were made to Turner's dynamic analysis [2], the Conel technique [1], and a modified version of the Conel technique that Spang and Chapple developed [7]. The overall result was that all three methods produce shortening and extension axes with roughly the same orientation, and the strain is substantially underestimated where magnitudes were available from Conel and modified Conel techniques.
Application of the strain gauge technique to field data produced new questions: First, how wide is a twin set? Under the microscope, thick twins have a measurable thickness of twinned material. Thin twins (formerly known as untwinned lamellae or microtwins), on the other hand, resemble dark lines crossing the grain. The field samples I was working with had many thin twins and a few thick twins. Conel [1] had conducted a thorough analysis of the optics and concluded that the visible line was the result of an extremely thin twin. Therefore, the line represents a twin that is presumably equal to or thinner than the width of the line. Thick twins are bounded by lines of finite width, creating the question of whether the line width should be included as part of the twin thickness. The other question is whether or not the twin sets actually record the correct amount of strain. If twin sets are passive markers, each one should have exactly the correct strain for its orientation. However, at grain boundaries, where a twin set hits the adjacent grain, it might be prevented from deforming enough, causing the twin to underestimate the bulk strain. There might be other reasons that the twins under-record the strain. Neither the width nor accuracy questions can be definitively resolved using field data, but it was necessary to try.
In my field data, the ratio of thin twin-line thickness to actual twin thickness was the best fit by a ratio R = 0.5. A ratio of 0.75 made many thin twins thicker than thick twins, and a ratio of 0.25 generated too many very small-strain twin sets when compared to the simulated data. The inner thickness of the thick twins was used based on the premise that the black border was entirely an optical effect due to the abrupt change of crystal orientation. The accuracy of the principal strain magnitudes was guesstimated by comparing the strain calculated from twinning in a very gently folded bed close to an unfolded bed containing deformed brachiopods. The twin strain in the fold was taken as the average of the inner arc and outer arc values in order to average out the bending strain. The boundary-value strain for the fold is then the curved-bed shortening plus the internal strain. This latter concept was tested on a finite element fold model from Dieterich [8] and matched to within 17%. The final strain for the outcrop fold was −3.6%, the brachiopod strain was −7.6%, and the shortening directions were 30 • from one another. I concluded that the magnitude accuracy of the strain gauge should be no worse than a factor of 2. It is obvious that a large number of assumptions have been required to arrive at this point, leaving the accuracy of the strain gauge open to question. Consequently, I held off publishing my field results until 1975 [9], after the experimental confirmation of the accuracy of the technique, discussed next.

Experimental Confirmation
I finished my dissertation in 1971 and published the strain gauge technique in 1972, but being unsure about the accuracy of the method when applied to real data, I held off publishing my results for the natural folds I had studied. After hearing a presentation of my results at a meeting, Mel Friedman and John Logan of the Texas A&M University Center for Tectonophysics group offered to let me try the method on thin sections from a series of experiments that John had run. In the experiments, cylindrical samples of limestone were shortened to axial strains of 0 to −8.4% (shortening negative). Then teaching at Syracuse University, I wrote and received a National Science Foundation grant to measure the samples and compare the strain gauge results to the experimental values. Later, after I became acquainted with John Logan personally, he told me that he and Mel had reviewed my grant proposal for the NSF and did not think the method would work but believed that I should be allowed to try it. This is an admirable mindset for grant reviewers, for which I thank them. The grant also helped to support an undergraduate research project by Larry Teufel at Syracuse in which he measured twin strain in a naturally deformed rock. Larry later studied at Texas A&M, where he ran additional experiments to further develop the strain gauge technique, as will be described later.

Methods
The experimentally deformed rock was the fossiliferous and spar-cemented Indiana limestone [10]. Cylindrical samples were deformed at 100 MPa confining pressure, room temperature, and at a constant strain rate of 10 −4 per second. The actual porosity of the individual samples before deformation is unknown. On average, the undeformed limestone had a porosity of about 13%, which was reduced by the deformation to between 5.1% and 0.4%, depending on the sample. The load history began with the application of confining pressure, followed by axial shortening. Two samples were subjected only to confining pressure, which was applied to the side of the jacketed cylinder while zero displacement was maintained by the axial piston. This induced porosity loss, which caused twins to form. These samples provide insight into the effect of two-dimensional compaction. Six additional samples were first loaded laterally as described and then shortened using the axial piston.
As a side note, while preparing the current paper, I realized that the captions on Figures 1 and 2 in my 1974 paper [10] are reversed. That is confusing, even for me. As another side note, in my original papers, I spelled gauge as gage because that how is my first text on solid mechanics spelled it [11] but changed to gauge when I saw that other geologists and rock mechanics generally spelled it that way.
The methods tested using the experimental data are the Groshong [5] strain gauge technique, the original Conel [1] technique, the Spang and Chapple [7] modified Conel technique, and the then-new Spang [12] numerical dynamic analysis method. All except the last method depend on the thickness of the twins. The Spang [12] numerical dynamic analysis (NDA) method assigns a shear magnitude of 1 to every C and T axis, uses the tensor transformation law to transform all the data into thin-section coordinates, as did Conel [1], and then averages the components. The eigenvalues and eigenvectors of this average tensor give the C (NDA 1 ) and T (NDA 3 ) axes plus the intermediate axis (NDA 2 ).
Geosciences 2021, 11, 296 5 of 12 One advantage of this method over contouring is that it forces the C and T axes to be orthogonal, which contouring of the axes on a stereonet does not necessarily accomplish.
ure 2). The values that were selected, based primarily on the calculated statistical error, were the ratio R between measured and actual thin twin width of 0.5 and the inner width of the thick twins. In hindsight, today I might choose the outer thick twin width because it provides a better fit to the experimental strain magnitudes (Figure 2). Note that the maximum strain magnitude difference is about 2. It seems unlikely, even today, that the difference between, for instance, 6.5% and 8.5% shortening would significantly impact any interpretive hypothesis. In a later series of experiments Teufel [13] found that the strain gauge calculation consistently slightly underestimated the experimental strain magnitude, another reason to slightly increase the thicknesses assigned to thick twins and perhaps even to the thin twins. Two representative samples will be discussed in detail; sample 282, which underwent only lateral compaction, and sample 278, which underwent a large axial shortening in addition to lateral compaction. Sample 282 gives the clearest picture of the effect of twodimensional compaction because its porosity was reduced to almost zero (0.8%) by the initial loading. Even though the axial displacement is zero, this form of loading produces a deviatoric extension in the axial direction because of the lateral shortening due to porosity loss. The strain gauge technique indicates a deviatoric extension with an axis about 10° from the cylinder axis (Figure 3a), a fairly accurate estimation of the direction. The deviatoric extension magnitude is +1.3%. This is comparable to the experimental value, which could be as much as +4.3% depending on the exact initial porosity, which is unknown. The Turner [2] dynamic analysis extension axis distribution (Figure 3a) appears to be random, and the compression axes are equally random [10]. Even though grains are clearly deformed by twinning, there is no obvious directional information given by the Turner Figure 2. Computed vs. experimental maximum shortening strain (ε 1 ) as a function of different thicktwin width measurements, all data (modified from [5]). Sample numbers (vertical) are followed by the experimental axial shortening in percent. A dashed line represents perfect agreement. Samples 280 and 282 were maintained at zero axial shortening while being loaded by confining pressure.

Results
Due to the finite width of the twin boundary lines, alternative measures of twin width were again tested to find the values that best fit the experimental strain magnitudes ( Figure 2). The values that were selected, based primarily on the calculated statistical error, were the ratio R between measured and actual thin twin width of 0.5 and the inner width of the thick twins. In hindsight, today I might choose the outer thick twin width because it provides a better fit to the experimental strain magnitudes (Figure 2). Note that the maximum strain magnitude difference is about 2. It seems unlikely, even today, that the difference between, for instance, 6.5% and 8.5% shortening would significantly impact any interpretive hypothesis. In a later series of experiments Teufel [13] found that the strain gauge calculation consistently slightly underestimated the experimental strain magnitude, another reason to slightly increase the thicknesses assigned to thick twins and perhaps even to the thin twins.
Two representative samples will be discussed in detail; sample 282, which underwent only lateral compaction, and sample 278, which underwent a large axial shortening in addition to lateral compaction. Sample 282 gives the clearest picture of the effect of two-dimensional compaction because its porosity was reduced to almost zero (0.8%) by the initial loading. Even though the axial displacement is zero, this form of loading produces a deviatoric extension in the axial direction because of the lateral shortening due to porosity loss. The strain gauge technique indicates a deviatoric extension with an axis about 10 • from the cylinder axis (Figure 3a), a fairly accurate estimation of the direction. The deviatoric extension magnitude is +1.3%. This is comparable to the experimental value, which could be as much as +4.3% depending on the exact initial porosity, which is unknown. The Turner [2] dynamic analysis extension axis distribution (Figure 3a) appears to be random, and the compression axes are equally random [10]. Even though grains are clearly deformed by twinning, there is no obvious directional information given by the Turner technique. The numerical dynamic analysis (NDA) locates a deviatoric extension axis (NDA 3 ) within about 20 • of the cylinder axis, a reasonable approximation, given that the T axes appear to be random. technique. The numerical dynamic analysis (NDA) locates a deviatoric extension axis (NDA3) within about 20° of the cylinder axis, a reasonable approximation, given that the T axes appear to be random. Sample 278 is one of the most deformed, having an axial shortening of −7.5%. The computed strain gauge maximum shortening direction is within about 2° of the experimental direction, and the computed strain magnitude is just 0.5% larger than the experimental value, a very good match. For this sample, Turner compression axes show a clear clustering around the cylinder axis ( Figure 3b) and the extension axes form a girdle perpendicular to the cylinder axis (see [10]). NDA1 (average C axis) is also almost exactly parallel to the cylinder axis. The NDA1 magnitudes are 0.5% or less ( Figure 4a) and do not correlate with either the stress or strain magnitudes.
The twinning strain gauge does not measure volume change, therefore the calculated strains represent the deviatoric component of the total. For deviatoric strain, two principal strain magnitudes must have one sign and the third will have the opposite sign because the sum of the three must be zero. The odd-signed axis is always the best located in the experiments because of the cylindrical symmetry of the strain; it is in the direction of the cylinder axis. In the experiments, the same-signed axes should be identical in magnitude, and therefore their directions are not well constrained.

Noise Reduction
Noise is defined as (1) the size of the statistical standard error [5], (2) a negative expected value for the strain of a twin set as calculated from the best-fit result [10], and (3) the size of the difference between the expected value for each twin set and the twin strain, Figure 3. Comparison between sample 282 with zero axial shortening and sample 287 with −7.5% axial shortening. Turner [1] dynamic analysis axes (circles and + signs), along with calculated principal axes for the Spang [12] numerical dynamic analysis (NDA) and the Groshong strain gauge technique (ε); all data modified from [10]. Lower-hemisphere equal-area projections in the thin-section coordinate system. The orientation of the cylinder axis relative to the thin section coordinate system is shown schematically by rectangles. The long axis of the sample cylinders (short lines labeled Cyl. Axis) are not quite parallel to the long axis of the thin section: (a) sample 282; (b) sample 278.
Sample 278 is one of the most deformed, having an axial shortening of −7.5%. The computed strain gauge maximum shortening direction is within about 2 • of the experimental direction, and the computed strain magnitude is just 0.5% larger than the experimental value, a very good match. For this sample, Turner compression axes show a clear clustering around the cylinder axis ( Figure 3b) and the extension axes form a girdle perpendicular to the cylinder axis (see [10]). NDA 1 (average C axis) is also almost exactly parallel to the cylinder axis. The NDA 1 magnitudes are 0.5% or less ( Figure 4a) and do not correlate with either the stress or strain magnitudes.
The twinning strain gauge does not measure volume change, therefore the calculated strains represent the deviatoric component of the total. For deviatoric strain, two principal strain magnitudes must have one sign and the third will have the opposite sign because the sum of the three must be zero. The odd-signed axis is always the best located in the experiments because of the cylindrical symmetry of the strain; it is in the direction of the cylinder axis. In the experiments, the same-signed axes should be identical in magnitude, and therefore their directions are not well constrained.

Noise Reduction
Noise is defined as (1) the size of the statistical standard error [5], (2) a negative expected value for the strain of a twin set as calculated from the best-fit result [10], and (3) the size of the difference between the expected value for each twin set and the twin strain, as calculated from the best-fit result [10]. Potential sources of noise can include measurement error, inhomogeneous strain, or a complex strain history. My 1974 experimental results were considerably "noisier" than the results from my field data, a fact that needed an explanation.
Can the sources of the noise be identified and then reduced using methods that will also work with field data? In 1974 [10], I inferred that much of the experimental noise probably resulted from inhomogeneous strain due to porosity collapse. This was later confirmed experimentally by Teufel [13] who, using the same material (Indiana limestone), obtained much less noisy results by not including twin sets associated with porosity collapse. These grains were identified as ones where the twins start at grain boundaries but do not continue across the whole grain. This criterion had been put forward by Friedman [14] and was applied by Teufel [13]. Grains similar to this were included in the data sets in my 1974 paper, none of which was recorded, making later removal impossible.
Removing the largest deviations with respect to the expected values worked well with my 1974 experimental data. This cleaning technique should remove significant measurement errors and reduce the effect of an inhomogeneous strain. Removing the 10% largest deviations seemed to provide a reasonable balance between removing inconsistent measurements but without eliminating too much of the hard-won data [10]. This strategy also worked well for experiments by Teufel [13] and Friedman et al. [15].
Negative expected values pose a different kind of problem. Negative expected values indicate the wrong sense of shear for the orientation, indicating a grain that should not be twinned given the calculated strain tensor. Results in my 1974 paper [10] show that some negative expected values can arise from inhomogeneous strain. It seems unlikely that the strain would be perfectly homogeneous throughout every grain in every direction. Each grain has only three directions in which it can deform by twinning, and adjacent grains must provide some resistance to the formation of twin sets that try to deform a boundary that is not suitably oriented to twin. Thus, some twin sets with small negative expected values are no surprise, even in rocks with no initial porosity. However, large numbers of negative expected values would seem to imply a different process.

Multiple Deformations
The experimental study by Larry Teufel [13] addressed the question of multiple deformations as a cause of negative expected values. In his experiments, for single and two coaxial deformations, the percent negative expected values were 4 to 14%, the calculated strain magnitudes were within 10% of the experimental values, and the odd-signed axis was within 3 • of the experimental value. He found that obliquely superimposed strains were indicated if the calculated results had more than 40% negative expected values. When separately calculating the strain tensor for the positive and negative expected values, if the separated results have lower standard errors than for the combined data, then two deformations with the odd-signed strain axes approximately orthogonal to one another is indicated. The errors in magnitude and direction of the principal strains are about 20% and 6 • , respectively. When the superposed strains are at 45 • to one another, and positive and negative expected values are calculated separately, the errors in the two calculated strain tensors are about four times as large as for the single and coaxial strain experiments. The ability to calculate and interpret negative expected values turns out to be one of the significant benefits of using the strain gauge technique.

How Much Data Are Needed?
Measuring twins on a U-stage is a slow process, and some people may find it to be a bit tedious. Therefore, an important practical question is how much data are required for a trustworthy result. On the U-stage, one thin section does not access all possible twin orientations; two perpendicular sections should make it possible to observe almost all twin directions, and three would guarantee the observation of every possible direction, but are they all necessary? Additionally, how many twin sets must be measured for a good strain result? In 1984, we [16] addressed these questions using two outcrop samples and one experimental sample. We measured 150 or more grains total from three perpendicular thin sections and then calculated the strain from subsets of various sizes. The data were cleaned by removing the 20% largest deviations. We concluded that 25 twin sets from each of 2 perpendicular thin sections were sufficient. With these data, the magnitude of the odd-signed axis should be within 10% of the true value and within ±6 • −8 • in orientation.
Using just 25 twin sets from one thin section yielded a precision of about 20% in strain magnitude and ±10 • to 12 • in the angle. This should be satisfactory for most purposes.

Stress or Strain?
There has long been a tendency to refer to dynamic analysis [2] and related techniques as stress analysis, whereas methods using measured twin widths are referred to as strain analysis. I suggest that, if you can visually detect it, it is strain. For natural deformation, stress must be calculated from the strain, and that requires knowing the appropriate flow law and its physical constants. As shown previously, almost any reasonable calculation or contouring technique will provide a good estimate of the odd-signed strain axis direction, equivalent to the corresponding stress axis if the strain is irrotational. Two different approaches for obtaining stress magnitude from twinned calcite are discussed next.
The standard approach to finding the stress magnitude from strain is via the appropriate flow law. The stress-strain relationship for the experimentally deformed samples from 1974 is given in Figure 4a. The relationship is probably nonlinear, but a straight-line best fit is shown. The NDA 1 values are not correlated to the strains but are most likely related to the strength of the C and T axis concentrations. The twin strain calculated by the strain gauge technique is well correlated to the stress except for the largest strain value. Here, the twin data have had the 10 largest deviating twin sets removed from the calculation. This improves the agreement between the experimental and calculated values by a small amount except for the maximum strain magnitude, which remains underestimated. This is another reason why it probably would have been better to use the outer thick twin width in the calculation (Figure 1). The stress-strain relationship in Figure 4a is not predictive for field use unless the flow law is identical to that in the experiments.
but are they all necessary? Additionally, how many twin sets must be measured for a good strain result? In 1984, we [16] addressed these questions using two outcrop samples and one experimental sample. We measured 150 or more grains total from three perpendicular thin sections and then calculated the strain from subsets of various sizes. The data were cleaned by removing the 20% largest deviations. We concluded that 25 twin sets from each of 2 perpendicular thin sections were sufficient. With these data, the magnitude of the odd-signed axis should be within 10% of the true value and within ±6°−8° in orientation. Using just 25 twin sets from one thin section yielded a precision of about 20% in strain magnitude and ±10° to 12°in the angle. This should be satisfactory for most purposes.

Stress or Strain?
There has long been a tendency to refer to dynamic analysis [2] and related techniques as stress analysis, whereas methods using measured twin widths are referred to as strain analysis. I suggest that, if you can visually detect it, it is strain. For natural deformation, stress must be calculated from the strain, and that requires knowing the appropriate flow law and its physical constants. As shown previously, almost any reasonable calculation or contouring technique will provide a good estimate of the odd-signed strain axis direction, equivalent to the corresponding stress axis if the strain is irrotational. Two different approaches for obtaining stress magnitude from twinned calcite are discussed next.
The standard approach to finding the stress magnitude from strain is via the appropriate flow law. The stress-strain relationship for the experimentally deformed samples from 1974 is given in Figure 4a. The relationship is probably nonlinear, but a straight-line best fit is shown. The NDA1 values are not correlated to the strains but are most likely related to the strength of the C and T axis concentrations. The twin strain calculated by the strain gauge technique is well correlated to the stress except for the largest strain value. Here, the twin data have had the 10 largest deviating twin sets removed from the calculation. This improves the agreement between the experimental and calculated values by a small amount except for the maximum strain magnitude, which remains underestimated. This is another reason why it probably would have been better to use the outer thick twin width in the calculation (Figure 1). The stress-strain relationship in Figure 4a is not predictive for field use unless the flow law is identical to that in the experiments. Points are labeled with sample numbers: (a) differential stress (∆σ) vs. strain (ε) curve for Indiana limestone deformed at 100 MPa confining pressure, room temperature, and a constant strain rate of 10 −4 /s (compiled from data in [10,17]. The dashed line is the least-squares best fit to the experimental stressstrain points. The twin strain data have had the 10 largest deviating twin sets removed before calculation. Results from one sample are aligned vertically; (b) stress estimated from percentages of twinned grains, compiled from data in [17]. The twin-set percentage curves (blue) are cumulative; the one twin-set curve includes all grains with two and three twin sets. Experimental stress magnitudes are red squares and, in 276 and 278, identical.
After finishing my 1974 paper, I passed the thin sections along to Jamie Jamison and John Spang at the University of Calgary. As part of his MSc thesis, Jamison developed a very different approach to estimating stress magnitude. The basic concept is that, because there is a critical resolved shear stress required to initiate twinning, in a randomly oriented crystal aggregate the greater the differential stress applied to the aggregate, the greater the number of twin orientations that will be activated. Jamison and Spang [17] derived a relationship between the expected percentage of twinned grains and what they termed the resolved shear stress coefficient, S 1 , which is a function of the angles between the pole to the twin plane, the glide direction, and the principal stress directions. Figure 4b is a plot of the percentage of grains having one, two, and three twin sets. These percentages are measured on a thin section and plotted, giving the value of S 1 for that sample. S 1 is converted to stress with the equation ∆σ = t c /S 1 , where t c is the critical resolved shear stress for twinning. They used a value for t c of 10 MPa based on the research by Turner, Griggs, and Heard [18], and Friedman [19].
I have plotted the Jamison and Spang measurements [17] for their samples in Figure 4b. For an internally consistent result, the cumulative percent of one, two, and three twin sets for a thin section should fall on a vertical line. The relationship between ∆σ and S 1 is strongly nonlinear, with ∆σ going to infinity as S 1 goes to zero. Thus, the resolution of ∆σ using the graph as presented is not very good for S 1 < 0.1, equivalent to ∆σ > 100 MPa. For the two least deformed samples (264, 262), the lines are close to vertical and the predicted stresses are very close to the experimental values. The other samples provide stress estimates that are too high. The worst fit is the most deformed sample (261) for which the line between two and three twin sets is the farthest from vertical and the stress estimate is significantly greater than the experimental value (the strain gauge result for this sample (Figure 4a) is the biggest underestimate of the strain magnitude). The curve for three twin sets (Figure 4b) yields the most accurate result, as noted by [17]. Using the three twin-set values, the predicted maximum compressive stress in every sample agrees with the experimental value to within 30 MPa.
The stress-strain relationships in Figure 4a,b may be misleading as estimates of the stresses in naturally deformed rocks because the results depend on the load history. The experiment was conducted at a constant strain rate, which causes the stress to increase with time under the given conditions of temperature and pressure. As noted by [17], Friedman and Heard [20] found that the number of twins and the number of twinned grains increases with time at constant differential stress in a creep (constant stress) test. The effect of the loading rate on the stress and twin strain, along with the temperature and pressure effects, should be investigated carefully before any attempt is made to infer stress magnitudes from twin strain.

Discussion
It is now possible to discuss the interpretation of a representative fold from the group of natural folds that I originally studied [9]. The example in Figure 5 is selected because it is the best geometric match to the finite-element fold model of Dieterich [8]. The values of shortening at corresponding locations in the natural and model folds differ by an order of magnitude or more. My initial conundrum was whether or not the calculated strain could be incorrect by such a large amount. The experimental results show that the strain gauge calculation is accurate to 20 or 30% of the magnitude value; thus, the order-of-magnitude discrepancy is real. This is also consistent with the fabric of the limestone itself because the fossils are not obviously deformed, indicating that the maximum shortening at the grain scale is almost certainly no more than 10%.
Two additional deformation mechanisms are obvious in the fold ( Figure 5): tension fractures, now filled with secondary calcite, and wavy stylolites containing insoluble residue of quartz, dolomite and, in some, presumed organic material. Stylolites are abundant throughout most of the layer but cannot be positively identified near the top because there the layer consists of fine-grained dolomite without fossils. The number of stylolites shown is the minimum because zones of insoluble residue that do not clearly truncate a fossil are not included. The stylolites are roughly normal to the measured shortening directions and to the predicted shortening directions in the model. Some stylolites curve toward bedding in the outer arc in both the anticline and in the right-hand syncline, as predicted by the model. The filled tension fractures at the top and base of the layer represent outer-arc extension. The large fracture on the right limb of the anticline initiated at its lower, most nearly horizontal boundary and opened obliquely with a left-lateral sense of shear, as indicated by the fabric of the calcite filling [9]. The fracture orientations are parallel to the predicted maximum compressive stress directions in the buckle fold model [22]. All three deformation mechanisms are contemporaneous. The strain is partitioned among them according to their individual flow laws and each has a different effect on the final geometry. Presumably, only a small fraction of the total strain is partitioned into calcite twinning with most of the strain being by pressure solution. Duplicating the results in Figure 5 is a worthy challenge for anyone conducting mechanical modeling.
Geosciences 2021, 11, x FOR PEER REVIEW 10 of 12 discrepancy is real. This is also consistent with the fabric of the limestone itself because the fossils are not obviously deformed, indicating that the maximum shortening at the grain scale is almost certainly no more than 10%. Figure 5. Comparison between calcite twin strains in a portion of the medium-dip fold layer, modified from [9] and the inset finite-element buckle fold model [8]. Both show shortening strains (negative) in percent. In the field example, arrows are maximum shortening directions in the plane of section, shaded areas are extension fractures filled with secondary calcite, and pale gray lines are stylolites. Every stylolite shown truncates a fossil somewhere along its length. In the inset model [8,9], short lines are principal shortening directions with magnitudes provided by Jim Dieterich [21]. Strains in the model are about an order of magnitude larger than in the field example.
Two additional deformation mechanisms are obvious in the fold ( Figure 5): tension fractures, now filled with secondary calcite, and wavy stylolites containing insoluble residue of quartz, dolomite and, in some, presumed organic material. Stylolites are abundant throughout most of the layer but cannot be positively identified near the top because there the layer consists of fine-grained dolomite without fossils. The number of stylolites shown is the minimum because zones of insoluble residue that do not clearly truncate a fossil are not included. The stylolites are roughly normal to the measured shortening directions and to the predicted shortening directions in the model. Some stylolites curve toward bedding in the outer arc in both the anticline and in the right-hand syncline, as predicted by the model. The filled tension fractures at the top and base of the layer represent outer-arc extension. The large fracture on the right limb of the anticline initiated at its lower, most nearly horizontal boundary and opened obliquely with a left-lateral sense of shear, as indicated by the fabric of the calcite filling [9]. The fracture orientations are parallel to the predicted maximum compressive stress directions in the buckle fold model [22]. All three deformation mechanisms are contemporaneous. The strain is partitioned among them according to their individual flow laws and each has a different effect on the final geometry. Presumably, only a small fraction of the total strain is partitioned into calcite twinning with most of the strain being by pressure solution. Duplicating the results in Figure 5 is a Figure 5. Comparison between calcite twin strains in a portion of the medium-dip fold layer, modified from [9] and the inset finite-element buckle fold model [8]. Both show shortening strains (negative) in percent. In the field example, arrows are maximum shortening directions in the plane of section, shaded areas are extension fractures filled with secondary calcite, and pale gray lines are stylolites. Every stylolite shown truncates a fossil somewhere along its length. In the inset model [8,9], short lines are principal shortening directions with magnitudes provided by Jim Dieterich [21]. Strains in the model are about an order of magnitude larger than in the field example.
New rock deformation tests could answer many questions. Each new experiment will provide the opportunity to use every calculation technique now available to observe how they compare. A high priority should be to find out whether constant-stress (creep) experiments produce the same twin strains as constant strain rate experiments. The results should show whether or not stress magnitudes can be predicted from the twin strain. It would also be of interest to see if the final fabrics differ between creep and constant strain rate experiments at the same total strain. All experiments should be run in both wet and dry conditions to investigate whether water affects the strain partitioning.
Natural and human-made compaction is a significant issue in reservoir development for areas ranging from hydrocarbon production to waste disposal, to potable water production. Calcite twin strains might provide new insights into the process. Experiments should begin with uniaxial compaction tests, representing compaction due to the weight of the overburden. Is it possible to distinguish between natural compaction and human-made compaction? Twins that do not cross an entire grain are thought to indicate the collapse of adjacent pores. What would happen if the strain for only such grains is calculated?
There are a variety of other interesting questions that could be addressed experimentally. For example, how critical is the magnitude of the difference between the principal strain magnitudes? This could be examined with true triaxial experiments in which the three principal strains can be controlled independently. It would also be interesting to see the results of the deformation of samples with a strong initial crystal-orientation fabric. Can the imposed strain be accurately measured or does the rock fabric prevent it? I hope someone finds one or more of these suggestions of interest and takes action.

Conclusions
The experimental results show that the odd-signed principal axis orientation is the most accurately located. In cases in which one principal strain magnitude is appreciably different from the other two, any calculation technique will be effective in finding that axis, generally to within about 10 • [10]. The strain gauge technique gives the correct strain magnitude to about ±20-30% of the magnitude. Experiments by [15] later confirmed the accuracy of 1% strain over the range of 1-12% axial shortening. Experiments by [13] showed that non-coaxial deformation can be recognized if the percentage of negative expected values of strain in the grains is greater than 40%. Calculating the positive and negative expected values separately can give the correct superposed strain magnitudes and directions. The stress analysis technique of [17] predicts the experimental stress magnitude to better than 30 MPa using their resolved shear stress coefficient vs. percentage of twinned grains relationship but only for grains with three twin sets. The stress axes will coincide with the strain axes only if the strain is irrotational.
Funding: This research received no external funding.