A Mathematical Technique for Analyzing Folds With the Computer Program “FOLDPI”

A mathematical technique for analyzing folds was proposed instead of the tedious and slow graphical method. Procedure of this technique comprises converting the data of bedding planes to pole attitudes, calculation of the mean pole vector of fold limbs, obtaining the best fit π-circle, determining the fold geometric properties and finding fold cylindricity. This procedure was carried out by FOLDPI, a GWBASIC computer program written for the purpose of this application. Most of the geometrical properties of fold were dealt with. In addition, an example taken from the Sinjar Anticline was used for testing the validity of this technique. The results of testing the program against manually obtained solutions proved that this technique can be very helpful in getting faster and more accurate results ــــــــــــــــــــــــــــــــــــــــــــــــــ يبوساحلا جمانربلا مادختساب تايطلا ليلحتل ةيضاير ةقيرـط FOLDPI


INTRODUCTION
Conventionally, folds were analyzed using the common graphical technique of πdiagram, which is one of the stereographic projection applications. This graphical method has a worldwide usage and it has an advantage of graphically showing fold geometry but it is tedious and time consuming specially in plotting and counting the S-poles on the stereonet. Recently, and for the sake of faster and easier techniques, many structural geologists have attempted to modify their related methods towards the trend of mathematics and computer programming approaches. Accordingly, the present author suggests a mathematical technique for digital execution of this π-diagram and the determinion of the geometric parameters of folds using mathematics and computer program.
Previously, some authors made contributions in this trend. They performed some steps in this respect. Ramsay (1967) suggested two mathematical techniques in the scope of fold analysis. The first was applied for determination of unimodal poles distribution by vectors of directional cosines. While the second method was used for determining the best-fit π-circles of cylindrical folds. Bengston (1980) mentioned, marginally, about this idea through out his study of tangent diagram. Ramsay and Huber (1987) described the methods of Ramsay (1967) by " the accuracy contained in such methods is only justified if it is imperative to exploit the full potential of very precise primary data".
π-diagram of fold analysis is a method applied by using an equal area stereographic projection (stereonet) to plot the perpendiculars to the bedding planes (S-poles). Attitudes of such bedding planes are collected along traverse which must be transverse to the fold axis. Plotting the S-poles of such a fold (point diagram), counting the concentration of points in 1% of the stereonet area, and lastly contouring the counted values in the form of contour percentages (contour %) to produce the π-diagram of this fold. The advantage of this diagram is the geometric view of the fold parameters. Such parameters are fold axis, fold plunge, fold symmetry, interlimb angle and attitude of axial plane.
In π-diagram, if the fold is perfectly cylindrical, the bedding S-poles fall, perfectly, along a great circle of the stereonet (π-circle). When the fold transforms geometrically to non-cylindrical types, the scattering of the poles around π-circle becomes more pronounced. The perpendicular to the π-circle is called π-axis that is coinciding with the fold axis. Consequently, if the fold is perfectly cylindrical the angle between each S-pole and fold axis is perfectly 90° (Fig. 1). Practically, the S-poles of any natural fold do not lie exactly on a certain great circle but fall in zone around this circle. Nevertheless, the human error in field measurements plays a role in the accuracy of these measurements; which amount to ±2 degree (Ramsay and Huber, 1987). However, the mean reason responsible for the scattering of S-poles is the absence of perfectly cylindrical folds in the field. Many of natural folds are of cylindrical, sub-cylindrical and non-cylindrical styles. In this respect Ramsay (1967) designed his method for determining the best-fit π-circles for perfectly cylindrical folds. This is because only in this type of folds the π-axis makes a right angle with each S-pole; and he built up his mathematics on this property. It must be mentioned that in Ramsay (1967) terminology, the term cylindrical fold is analogous to perfectly cylindrical one. In the more recent synthesis, Ramsay and Huber (1987) classified the folds into perfectly cylindrical, cylindrical, sub-cylindrical and noncylindrical types (Fig. 2). Previously, folds classified as cylindrical and non-cylindrical and some time cylindroid (Fleuty, 1964). In the classification of Ramsay and Huber(1987) the perfectly cylindrical fold has its S-poles lie perfectly on the π-circle.
While if more than 90% of the S-poles fall within an angle of ±10° from the constructed π-circle the fold should be termed cylindrical. But if more than 90% of the data lie within ±20° of the π -circle the fold is then called sub-cylindrical fold. Folds with more than 10% of their S-poles falling outside the limit of ±20° are termed noncylindrical (Fig. 2). In addition to the scope of the present work, the author extended the method of Ramsay (1967) from its application only on perfectly cylindrical to include all types of folds described by Ramsay and Huber (1987). This extension was based on the styles of poles distribution a round the best-fit π-circle.

METHODOLOGY
The proposed mathematical method comprises the following procedures: 1st-Converting the data of bedding plane attitudes of both fold limbs (strike direction/dip amount when strike was taken clockwise from dip direction) to pole attitude (dip direction/dip amount) and finally to their corresponding directional cosine vectors (α, β & γ). 2nd-Calculation of the mean vector of each fold limb and hinge area by unimodal poles distribution method (summing method). 3rd-Obtaining the best-fit π-circle for these two means of fold limbs, or three means (two limbs and a hinge area) 4th-Determination of fold parameters which reflect its geometry. 5th-Determination of the cylindricity of folds according to Ramsay and Huber (1987).

1st-Determination of directional cosines :
Mathematical equations needed for the determination of directional cosines from attitudes of bedding planes are derived in this work. So the angles α, β & γ are determined from strike directions and dip amounts. The following steps explain this procedure.
1. Conversion of bedding plane attitudes (strike direction /dip amount) to pole attitude (dip direction / dip amount). So, dr = 90 + ds , dp = 90 -dpl cdr = 90 -dr Where dr and dp are the angles of dip direction and dip amount of the pole, ds and dpl are of strike direction and dip amount of the bedding plane and cdr is the complimentary angle of dr. 2. Determination of the directional cosines. These can fall into four cases. Each case represents one of the upper four quarters of the Cartesian coordinate. It must be mentioned that in the four cases the angle γ is always equal to (90 -dp) with negative sign (Fig.4). The first case: if the pole of any plane falls in the first quarter (Fig.3).
Determination of the angle α: In Figure (3), suppose the line OA in the triangle OAB is unity. Then AB = sin dp and OB = cos dp.

Determination of the angle β:
In the triangle OAB, OBA is a right angle, suppose the hypotenuse OA is unity (Fig.  3). Therefore, OB = Cos dp and AB = Sin dp The angle BCO in the triangle BCO was drawn to be 90°. So, Sin dr = CB / Cos dp so CB = Sin dr . Cos dp Cos dr = OC / Cos dp so OC = Cos dr . Cos dp ABC is a right triangle, So according to the Pythagorean theorem (AC)2 = (CB)2 + (AB)2 Then, (AC)2 = Sin2 dr . Cos2 dp + Sin2 dp Lastly It must be noted that this equation is similar to that equation of Cosα, accept the angle cdr is replaced by the angle dr. In this case, when the angle α fall in the first quarter it has a positive sign whereas the angle β is negative (Fig. 4).
The third case: If the pole of any plane falls into the third quarter. Then, drnew = drold -180 cdr = 90 -dr Similarly,

And
In this quarter, α has negative sign. While β is positive (Fig. 4).

2nd-The unimodal pole distribution method:
The unimodal pole distribution was suggested and described by Ramsay (1967). It is also called method of summing vectors. Considering the pole as a unit vector, the method is used for determining the mean vector of poles of any geologic planes. This is done after determining the angles α, β and γ of each pole (unit vector) with respect to the coordinate axes x, y and z respectively (Fig. 5A). These angles, in the present work, are an output of the previous procedure (mentioned in 1st). The x, y, and z components of the side of the vector box (Fig. 5B) for each measurement are then calculated by Cosα, Cosβ and Cosγ respectively. These directional cosines are determined for all poles and sums of the vector components are calculated (Σ cosα, Σ cosβ and Σ cosγ). These sums give the dimentions of the x, y, and z components of the total vector sum and the diagonal of this box gives the strength of the total vector sum (TVS) which is equal to: ) Therefore, the mean vector direction with respect to x, y and z axes are given by: In the course of this work, data that were taken from a fold can be differentiated into two or three concentrations. If the folds are of chevronic or mostly chevronic style, two concentrations of pole distribution are found. That is because the hinge is angular and there is no hinge area then the two concentrations representing the two limbs. Whereas, in the box type folds or near this shape, three concentrations are found. Two of them for the limbs and the third represent the hinge area. It means that each fold has two or three mean vectors. Consequently, the users of this method (Unimodal poles distribution) must process each concentration of fold poles and determine the mean vector of each one. These mean vectors (with the angles α, β and γ) can be used for determining the best-fit Cos β=[Cos 2 cdr.Cos 2 dp+1-Sin 2 dp-Cos 2 dp. Sin 2 cdr]/[2.Cos cdr.Cos dp] π-circle of such a fold. Mean vectors can be plotted manually by setting the angles α, β and γ on the stereonet or they can be processed mathematically to determine the fold πcircle.

3rd-Determination of best fit π-circle of folds:
Fold π-circle, as it mentioned above, is a stereographic best-fitted great circle to bedding S-poles. The method for determining the best-fit π-circles was suggested by Ramsay (1967). In this method, π-axis was determined mathematically by finding the normal to each bedding S-poles of perfectly cylindrical fold. Then π-circle can be drown, stereographiclly, perpendicular to this normal. The method was based on the usual technique of minimizing the squares of the deviations of the observed bedding S-poles from this surface.
Ramasy (1967; pp.18-20) designed this method for number of poles that fall perfectly on the π-circle, so this method is constrained for perfectly cylindrical folds only.
The present author made a simple modification for wider range of applications including cylindrical and sub-cylindrical fold, which are dominant in the field. Therefore, it is modified by obtaining the mean vector of large number of poles for the two fold limbs or two-fold limbs and hinge area by the unimodal pole distribution. Then, determining of the π-circle best fitted to the two or the three concentrations. By this way of drawing the π-cirle, limits of ±10º and ±20º from π-cirle can be plotted and type of fold cylindricity can be found (As in E). According to this modification many types of folds can be identified.
The method was described in (Ramsay, 1967), and it is summarized here as the following: Suppose α, β and γ are angles between π-axis vector of any fold and the three coordinate axes X, Y and Z respectively . Also, suppose a, b and c are angles between any bedding S-pole vector and the same coordinate axes (Fig. 6).
By the values and signs of these angles the Cartesian coordinate position of π-axis was found. Fold π-circle can be drawn considering this π-axis normal to it.

4th-Determination of fold geometry:
Most of the important geometric properties of folds can be determined from this technique. Such properties are fold axis, fold plunge, interlimb angle, fold symmetry, attitude of axial plane and fold cylindricity.

Fold axis:
The angles α, β and γ which were calculated in (3rd page 6 ) are used for determining the attitude of fold axis (dip direction / dip amount). From figure (7), Also from the same figure, the triangle OAC has the right angle OCA and the hypotenuse OA suppose to be unity. dp = 90 -γ when dp is the dip amount of fold axis cos γ = ( 1+ A 2 + B 2 ) -1/2 cos α = A (1 + A 2 + B 2 ) -1/2 cos β = B ( 1+ A 2 + B 2 ) -1/2 AC = sin β and OC = cos β Also in the triangle, OBA is a right angle. OB = cos dp and AB = sin dp Whereas in the triangle ABC AC2 = AB2 + CB2 CB = (sin2 β -sin2 dp)1/2 In the triangle OBC, CB2 = OC2 + OB2 -2. OC. OB.Cos dr Sin2 β -sin2 dp = cos2 β . cos2 dp -2 .cos β .cos dp . cos dr 2. cos β . cos dp . cos dr = cos2 β + cos2dp -sin2 β + sin2 dp When dr is the dip direction of fold axis. The signs of the angles α and β serve as indicators to show in which quarter of Cartesian coordinate the fold axis was fall (Table  1). If the fold axis falls in the first quarter, dr remains without change. Whereas, if it falls in the second, third and forth quarter, then 90°, 180° and 270° are added to the angle dr for obtaining its correct dip direction).

Fold plunge:
Fold plunge depends upon dip amount of fold axis (dp). When (dp) is equal to zero, the fold is nonplunging. While if (dp) is more than zero, the fold becomes plunging. And increasing of dp angle means increasing of degree of plunging.

Interlimb angle:
According to this procedure, interlimb angle can be calculated by adding the absolute value of γ1 (of the first limb) to the absolute value of γ2 (of the second one). The angles γ1 and γ2 are always having negative sign (Fig.4) then it must be considered their absolute values when the interlimb angle was determined. This calculation must be done along a plane perpendicular to the fold axis. Then concentrations of the two limbs must be rotated until the fold axis becomes horizontal and determination of the interlimb angle was done along the N-S vertical plane (Fig. 8). dr = cos -1 [ (cos 2 β+cos 2 dp-sin 2 β +sin 2 dp)/(2*cos β* cosdp)] Attitude of fold axis = dip direction (dr) / dip amount (dp)

Fold Plunge = dp of fold axis
Where adpl1 and adpl2 are mean dips of the two limbs

Folds symmetry:
Fold symmetry can be found by comparing adpl1 with adpl2 or γlimb1 with γlimb2. Resultant, if dpl1 is equal to dpl2. Then the fold is symmetrical, otherwise it is asymmetrical with the vergency being towards the limb which has greater pole dip angle (dp).

Attitude of axial plane:
This attitude will be possible if the attitude of axial trace is known (Ramsay and Huber, 1987). While if the fold is of parallel type, axial plane can be obtained by joining the bisector of interlimb angle with the fold axis. This is because in the parallel fold the axial plane always bisects the interlimb angle. In this work the axial plane determination was restricted to parallel fold type because other types are more complicated to be processed mathematically.
Axial plane attitude can be represented by its dip amount Apdp and dip direction Apdr.
Apdr coincide with the dip direction of gentle limb.

5th-Determination of fold cylindericity:
Ramsay and Huber (1987) classified the folds into perfectly cylindrical, cylindrical, sub-cylindrical and non-cylindrical, which were described previously in the introduction (Fig. 2).

Attitude of axial plane = Apdp / Apdr
Different natural folds show different properties or different π-diagram models. And it is very complicated to put a solution for each model. So a simplification was made for putting a general procedure for all these models. This is to revolve and rotate the data (mentioned below) to make a general solution for all the cases and keeping the entire relative geometrical relationships constant with each other but not with the coordinate axes.

Revolving and rotating fold data
Many authors suggested various methods for rotation of oriented data. Saha (1987) designed a FORTRAN program for rotation of data by transformation matrix. Al-Azzawi and Al-Jumaily (2000) proposed a mathematical procedure for rotation of joint planes relative to bedding rotation by trigonometric method.
For the sake of applying this idea, two stages are used. Firstly, revolving of data, horizonally, until direction of π-axis coincides with east direction of stereonet. This is done by adding (90 -dr) when dr less than 90°, and subtracting (dr -90) when dr more than 90° to or from the angle dr.
Secondly, rotation of all data around Y coordinate axis ( N-S line in the stereonet) until π-axis becomes horizontal; it means rotation angle (R) = dp of π-axis (Fig. 8). When π-axis rotated by the angle R , π-circle became vertical plane. This is because πaxis is always perpendiocular to the plane containing bedding S-poles (π-circle). This rotation can be done by transformation matrix. Many authors such as Arfken (1970) and Saha (1987) described this method. And it is summarizes here by the followings: 1-Determination of L, M & N components for each S-pole vector before rotation. L = cos dp . sin dr , M = cos dp . cos dr and N = sin dp 2-Multiplication of the components L, M & N by the transformational matrix that is responsible for rotation of vectors anticlockwise around Y-axis and through an angle R.
The resultant are new L, M & N components ( after rotation)which are used to determine the direction dr and dip angle dp of a pole vector.
Mathematically, this fold classification could not be applied without revolving and rotating data (as it mentioned above). Rotating all data until π-axis becomes horizontal standardize all natural cases in into one form. Figure (9) shows the rotated state of Figure  (2) which responsible for this classification. Figure (9) can be used to plot the revolved and rotated bedding S-poles of any fold and to find the type of this fold according to its cylindricity.

If dr < 90 then dr = dr + (90 -dr ) clockwise rotation
Wheras, if dr > 90 then dr = dr -(dr -90) anticlockwise dp = sin -1 N new and dr = sin -1 (L new / cos dp) or dr = cos -1 (M new / cos dp) Ramsay and Huber (1987) suggested limits for difining fold cylindricity (Fig. 2). These limits are ±10 and ±20 around π-circle of perfect cylindrical fold. These limits were divided the stereonet into fields; each field represents one of fold types. Figure (9) showed these limits after revolution and rotation. Empirically, the present author derived mathematical equations to deal with these limits during the present technique. And using Lagrangian interpolation method mentioned in (Al-Azzawi, 2004) does this. For more explanation, the author exhibited curves shown in figure (10) representing these limits S-pole directions of stereonet normally range from 1° to 360°. So that, direction of each S-pole can be tested and determined in which field it fall. Counting of these poles and obtaining the percentage of each field was done. Then the fold can be classifying as in the followings: 1-If all S-poles fall on the plane of perfectly cylindrical (π-circle) (Fig. 11), then it is called perfectly cylindrical fold. 2-If 90% of S-poles fall between plane-10 and plane+10 around the mean π-circle ( Fig.9) or on and above the standard curve of cylindrical fold (Fig. 10), then the fold named cylindrical type.
3-If 90% of them fall between plane-20 and plane+20 (Fig.9) or on and above the standard curve of sub-cylindrical fold (Fig.10), the fold classified as sub-cylindrical type. 4-Otherwise or when more than 10% of the S-poles lie outside plane-20 and plane+20 ( Fig.9) or fall below the standard curve of sub-cylindrical fold (Fig.10), the fold becomes non-cylindrical type. So, mathematically or by computer programming, the number of bedding S-poles that fall between these limits can be determined.
The procedure for analyzing fold was designed in GWBASIC computer program called FOLDPI (see the Appendix)

Tested sample:
Data for checking the validity of this technique has been taken from Sinjar Anticline. Al-Azzawi (1982) studied this anticline through three traverses. They named Gaulat, Sinjar and Jeribi. The first one was used for this test. The anticline in this traverse is, stratigraphically, comprise Shiranish formation of Upper Cretaceous age (the older formation). And it is followed by Aliji Fn. Consequently, the geometric properties of this fold are shown in a π-diagram of Figure  (11). This Figure showed that the mean pole attitudes of first and second limb are 340/74 and 182/43, attitude of fold axis is 087/5, it is asymmetrical fold with N vergency, plunging fold with 7 degrees, the interlimb angle is 120º, dip amount of axial plane is 76º toward SE and more than 90% of bedding S-poles are fall within ±10 around the mean πcircle then the fold classified as cylindrical fold. A computer program output for this analysis is exhibited in Figure (12) that shows not only encouraging result but it is more easy, fast and accurate than the graphical method specially when it done by computer scheme.