Computers & Geosciences 28 (2002) 467–472
QUADROFa program to estimate principal curvatures of folds$ Sait Ismail Ozkaya* Baker Atlas GOEScience, P.O. Box 15425, Manama, Bahrain Received 12 September 2000; received in revised form 8 July 2001; accepted 11 July 2001
Abstract Determination of curvature of structural surfaces is a method for estimating apertures and the fluid-flow potential of opening mode fractures that are associated with folding. Curvature is usually calculated within a sliding window over a structural surface and presented as a function of map location. The exceptions are where a pair of minimum and maximum curvatures is needed for an entire fold structure. QUADRO program calculates this by fitting a quadratic trend surface to the fold. The magnitude and direction of the principal curvatures are given by eigenvalues and eigenvectors of the second derivative for the quadratic surface. Once the minimum and maximum curvatures are calculated some additional structural indicators such as Gaussian curvature can be determined and the direction of fold-related fracture sets are identified. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Principal curvatures; Quadratic trend surface; Jacobian; Eigenvalues; Eigenvectors
1. Introduction Since Murray (1968) introduced the concept, curvature of structural surfaces has been used as a means of estimating apertures and predicting fluid-flow potential of fold-related opening mode fractures (e.g. Ericsson, 1998; Stewart and Podolski, 1998; Lisle, 1994). Curvatures are commonly calculated within a sliding window over a structural surface and presented as a function of map location. In some cases, however, a single pair of minimum and maximum curvatures for an anticline or a syncline may be needed to determine Gaussian curvature (Lisle, 2000) or direction of fold-related fracture sets (Stearns, 1968). The QUADRO program accomplishes this by fitting a quadratic trend surface to the fold. The magnitude and direction of the principal curvatures are calculated by finding $ Code available from server at http://www.iamg.org/ CGEditor/index.htm *Corresponding author. Tel.: +973-212-234; fax: +973-212345. E-mail address:
[email protected] (S.I. Ozkaya).
eigenvalues and eigenvectors of the second derivative for the quadratic surface. The choice of quadratic trend surface is related to the fact that the second derivative of a second-order surface is constant, which allows the calculation of a single set of principal curvatures for an entire fold structure. The principal curvatures and directions vary across a structural surface for higher order polynomials.
2. Procedure Depth (or elevation) measurements from many locations over a structural surface m are input in tabular form with columns as x (UTM East) y (UTM north) and z (depth from datum to surface). The data points may be on a regular grid or randomly distributed. The basic assumption is that the surface is continuous, and not disrupted by major faults. If major faults are present, each fault block would have to be analyzed separately but small faults may be ignored. A single anticline or syncline must be selected for analysis since a quadratic surface cannot be fitted into complex surfaces with multiple folds.
0098-3004/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 0 1 ) 0 0 0 8 5 - 1
468
S.I. Ozkaya / Computers & Geosciences 28 (2002) 467–472
The analysis starts with finding the equation of a second-order trend surface and correlation coefficient. Eigenvalues and eigenvectors of the second derivative for the trend surface are calculated next. The program does not carry out a significance test for the correlation coefficient. 2.1. Second-order trend surface A quadratic trend surface, f ðx; yÞ is fitted on the structural surface as described by Davis (1973). The general equation of a second-order trend surface is f ðx; yÞ ¼ a0 þ a1 x þ a2 y þ a3 x2 þ a4 xy þ a5 y2 :
ð1Þ
The coefficients of this equation are calculated by the well-known least square method, which strives to minimize the sum of squares of deviations of structural surface from the trend surface. The goodness-of-fit of the quadratic surface is calculated from (Davis, 1973) Sr R2 ¼ ; ð2Þ St where, Sr is the sum of squares due to regression and St is the total sum of squares. Sr is given by X 2 ðP YÞ #2 Sr ¼ ; ð3Þ Y# n St is given by P P X ð YÞ2 X ð YÞ2 : ð4Þ St ¼ Y2 n n Y# is the value of the quadratic surface for a given pair of x0 and y0 coordinates (Eq. (1)) Y# ¼ f ðx0 ; y0 Þ: ð5Þ Y is the actual value of the structural surface at the same sample location. The correlation coefficient between the quadratic trend surface and the actual fold is the square root of the goodness-of-fit i.e. pffiffiffiffiffiffi R ¼ R2 : ð6Þ The correlation coefficient provides a means to evaluate how well the quadratic trend surface fits the actual fold and how representative are the principal curvatures and directions. 2.2. Maximum and minimum second derivatives of a quadratic trend surface
derivative i.e. kDjf 00 j:
ð8Þ
Curvature is defined as a positive quantity. However, in order to differentiate upward convex and concave curves, we deviate from this definition and use the second derivative directly instead of its absolute value. The equivalent second derivative expression for a surface, f ðx; yÞ; is the matrix, M: 1 0 2 q f q2 f " # B @x2 qx qy C fxx fxy C B ð9Þ M¼ ¼B 2 C: @ q f fxy fyy q2 f A qx qy
qy2
The entries of the matrix, M, are second derivatives of the surface with respect to x and y coordinates. For a quadratic surface, the second derivatives are constant and hence M is independent of location. From Eq. (1) " # 2a3 a4 M¼ : ð10Þ a4 2a5 The maximum and minimum curvatures are given by the eigenvalues of this matrix, which are found by solving the following characteristic equation: Mk ¼ lk:
ð11Þ
Since M is a 2 2 matrix, its eigenvalues are found by solving the following quadratic equation:
f l fxy
xx ð12Þ
¼ 0:
fxy fyy l In terms of the coefficient of regression Eq. (1), this becomes:
2a l a4
3 ð13Þ
¼ 0:
a4 2a5 l The matrix, M, is symmetric, hence the solution yields two real values l1 and l2 : These are approximately equal to the maximum and minimum curvatures. Both l1 and l2 are positive for a doubly plunging anticline or syncline. Eigenvalues have opposite signs for a saddle. One of the eigenvalues is zero for a cylindrical fold. The non-zero eigenvalue is positive for an anticline and negative for a syncline. Both eigenvalues are equal with no distinct eigenvectors, when the curvature is spheroidal (a dome or basin-shaped structure). 2.3. Gaussian curvature
The curvature k of curve f is defined as k¼
jf 00 j ð1 þ
f 02 Þ3=2
:
ð7Þ
The bedding dip angle is usually small and the denominator of Eq. (7) is close to one. Therefore curvature can simply be represented by the second
The product of the two eigenvalues gives the Gaussian curvature, G: G ¼ l1 l2 :
ð14Þ
The Gaussian curvature of a fold indicates whether the surface has been subject to extension or not (Lisle,
469
S.I. Ozkaya / Computers & Geosciences 28 (2002) 467–472
2000). Cylindrical folds have not undergone surface lengthening and their Gaussian curvature is zero. A neutral surface can be defined for a fold with zero Gaussian curvature, whereas no such neutral surface exists for folds with non-zero Gaussian curvature. 2.4. Principal curvature directions The eigenvectors from the characteristic Eq. (13) point to the direction of maximum and minimum curvatures. The azimuth of maximum curvature y is given by k11 y ¼ arctan ; ð15Þ k12 where k11 and k12 are components of the eigenvectors corresponding to the largest eigenvalue. The direction of minimum curvature is orthogonal to the direction of maximum curvature. The principal curvature directions should yield orientation of types I, II and III fractures that are frequently associated with folds (Stearns, 1968). For a cylindrical fold, the axis is perpendicular to the maximum curvature.
in Figs. 1 and 2 for these two examples with magnitude and direction of curvatures summarized in Table 1. The directions of principal curvatures are also indicated on the trend surface maps. Correlation coefficients for the second-order trend surfaces are 0.87 and 0.95, respectively, indicating a high degree of reliability. Fig. 1a shows a dome-shaped double-plunging asymmetrical anticline. The Gaussian curvature for this fold is significantly greater than zero indicating that it is a non-cylindrical anticline. The second anticline (Fig. 2a) has a near-zero Gaussian curvature and is a cylindrical fold. 3.1. Estimation of fracture aperture and transmissivity Lisle (2000) stated that structural curvature reflects strain only when the Gaussian curvature is zero. For types II and III opening mode fractures (Stearns, 1968) that are associated with cylindrical folds, such as the second example above, it is possible to estimate the total fracture aperture per unit distance. From half thickness of a brittle layer, T; and radius of curvature, r, one can calculate total extension, at, as follows (Aguilera, 1995): T : 2r
3. Examples
at ¼
QUADRO was tested on two examples. Structural contour and quadratic trend surface maps are presented
The radius of curvature is nearly equal to the reciprocal of maximum curvature (largest
50
m
N
-13
0
35
m
N
-1
0m
0m
30
-1
30
-1
-1250 m
-1250 m
(a)
ð16Þ
5 km
(b)
5 km
Fig. 1. Structural contour map of asymmetrical dome (map a) and corresponding quadratic surface (map b). See Table 1 for magnitude and direction of curvatures, which are also shown on quadratic surface map (b). Gaussian curvature is large indicating that this structure is not cylindrical fold.
470
-55
00
.00
ft
S.I. Ozkaya / Computers & Geosciences 28 (2002) 467–472
N
N -6000.00 ft
-6000.00 ft
.00
00
-55 ft
ft
(b)
0.00
5 km
-600
ft
.00 ft
0.00
-6500
-600
(a)
5 km
Fig. 2. Structural contour map of plunging anticline (map a) and corresponding quadratic surface (map b). Table 1 gives magnitude and direction of curvatures, which are also indicated on map. Gaussian curvature is nearly zero, which shows that this particular anticlinal is cylindrical fold.
eigenvalue) i.e. 1 r¼ ; l1
ð17Þ
where l1 is the largest eigenvalue from Eq. (13). The average fracture aperture, aav ; is the total extension divided by fracture density, d; at ð18Þ aav ¼ : d Fracture density can be estimated from cores and borehole image logs. The maximum radius of curvature for the second example is 167 km. For a half thickness of 10 m, this yields a total aperture of 0.06 mm/m. For a fracture density of 2/m, for example, the average aperture is 0.03 mm. Based on Eqs. (18) and (19) one can also estimate fracture transmissivity from the maximum curvature of
a cylindrical fold. Fracture transmissivity is related to its aperture, a; as follows (Snow, 1969): T¼
a3 : 12
ð19Þ
For the fold in example 2, the fracture transmissivity is 54 mD m/m. Such aperture and permeability calculations, however, must be regarded as first-order approximations. As numerous factors can modify fracture aperture, e.g. cementation, so the relationship of curvature to fracture aperture (or density) has only limited validity. Curvature cannot be used to predict fracture attributes of shear slip fractures of monoclines (Cooke et al., 2000) or thrust ramp anticlines or faultbend folds (Salvini and Storti, 2001; Erickson et al., 2001) because flat limbs are subjected to a volumetric strain larger than the crest. Besides, highly fractured
S.I. Ozkaya / Computers & Geosciences 28 (2002) 467–472 Table 1 Principal curvatures and directions for example structures in Figs. 1 and 2 Example 1 Correlation coefficient Maximum curvature Corresponding radius Azimuth of maximum curvature
0.87 3.8 103 (1/km) 263 km 277 deg
Minimum curvature Corresponding radius Azimuth of minimum curvature Gaussian curvature
2.13 103 (1/km) 468 km 8.25 deg 8 106 (1/km2)
Example 2 Correlation coefficient Maximum curvature Corresponding radius Azimuth of maximum curvature
0.95 5.9 103 (1/km) 167 km 78.5 deg
Minimum curvature Corresponding radius Azimuth of minimum curvature Gaussian curvature
0.9 103 (1/km) 5125 km 169 deg 1.06 106 (1/km2)
flanks may be flat because unfolding follows folding as a thrust sheet moves over a ramp. Jamison (1997), in a study on the Monkshood anticline, showed that aperture is related to total volumetric strain, but not necessarily to curvature. Schultz-Ela and Yeh (1992) reported that curvature failed to indicate strain especially in gently folded rocks.
471
Acknowledgements Thanks are extended to both reviewers, in particular to Dr. Richard Lisle, for their careful evaluation of the manuscript and their remarks, which greatly helped to improve the paper.
Appendix A Running QUADRO QUADRO is a Unix Sun Sparc C program but Turbo-C on a PC platform can also compile it. The program expects input from two files: a parameter and a data file. The parameter file must be called ‘‘quadro. inp’’ and must include input file name, UTM coordinates of origin, the depth unit, east–west and north– south boundaries of the structure to analyze. An example of an input parameter file is shown in Table 2 Input data file must contain in tabular form, east and north UTM coordinates and structural elevation (positive downward). An example of an input file is shown in Table 3. In this example, the file has 5 columns. Columns 1 and 2 are coordinates and column 3 is the depth (m) to the structural surface. The QUADRO program outputs results both to a file and on screen. An example of an output file (filename lst) is shown in Table 4. Optionally, if the zcap flag is 1 (see Table 1) the program generates a trend surface file (with extension. zcap), which lists elevations to the trend surface at the same sample points as in the input data file. This file can be used to plot contour map of the
4. Discussion and conclusion Principal curvatures of second-order trend surfaces capture deformation characteristics of folds that are important for the analysis of fold-related fractures. QUADRO is designed to determine curvatures of an entire fold structure, which can be used to determine Gaussian curvature and to predict the aperture and direction of fold-related fractures. The program can be modified to calculate the direction and magnitude of curvatures within a sliding window as a function of location. The size of the window may be changed. A small window may be used to detect faults, whereas a large window can be used for predicting density variation of fold-related fractures. Fitting a quadratic surface before determining curvature helps to identify the significant structural features by eliminating small-scale variations and smoothing out localized complexities. Where the curvature is not constant and not well represented by a quadratic surface, it may be preferable to fit a higher-order polynomial surface to the structural fold.
Table 2 Example of parameter input filea input file name: Depth unit 0 m 1 ft: Number of columns in input file: x column: y column: z column: Missing value: Xleft: Ybottom: Scale: zcap file flag=1 generates x y zcap file: Xmin: Xmax: Ymin: Ymax: a
Shuaiba 0 5 1 2 3 9999 467 000 2 354 000 1000 1 467 000 477 000 2 354 000 2 365 000
If depth flag is 1, depth is converted to meters. If zcap flag is on, QUADRO program generates trend surface output file (see Table 5).
472
S.I. Ozkaya / Computers & Geosciences 28 (2002) 467–472
Table 3 Example of data input file. File must have tabular format. There must be one column for UTM East, one column for UTM North and one column for depth to structural surface. Column numbers for variables are specified in parameter input file UTM East UTM North Variable 1 Variable 1 Variable 2 467 000 467 100 ....... . 467 200 467 300 Endfile
2 365 000 2 365 000
1336.3 1334.0
18.3 18.1
0.21 0.42
2 365 000 2 365 000
1334.5 1334.3
18.2 18.0
0.63 0.34
Table 4 Example of data output file. QUADRO program outputs correlation coefficient of second-order trend surface, magnitude and direction of principal curvatures and corresponding radius of curvature Input parameters: File name X origin (UTM) Y origin (UTM) Scale
Shuaiba 467000.00 2354000.00 1000.00
Results: Goodness-of–fit Correlation coefficient Maximum curvature (1/km) Radius of maximum curvature (km) Azimuth of maximum curvature Minimum curvature (1/km) Radius of minimum curvature (km) Azimuth of minimum curvature
0.76 0.87 0.0038 263.47 277.11 0.0022 467.97 8.25
Table 5 Example of trend-surface output file. This file can be used to generate contour map of trend surface UTM East
UTM North
Depth (m)
467 100.00 467 400.00 ....... ....... . 467 700.00 468 000.00
2 364 900.00 2 364 900.00
1329.54 1324.27
2 364 900.00 2 364 900.00
1319.35 1314.76
trend surface. An example of a trend-surface file is shown in Table 5.
References Aguilera, R., 1995, Naturally Fractured Reservoirs. PennWell, Tulsa. 521pp. Cooke, M.L., Mollema, P.N.,Pollard, D.D., Aydin, A., 2000. Interlayer slip and joint localization in the East Kaibab monocline, Utah: field evidence and results from numerical modelling. In: Cosgrove, J. W., Ameen, M. S. (Eds.), Forced Folds and Fractures, Special Publication 169. Geological Society, London, pp. 23–49. Davis, J.C., 1973, Statistics and Data Analysis in Geology. Wiley, New York, 550pp. Erickson, S.G., Strayer, L.M., Suppe, J., 2001. Initiation and reactivation of faults during movement over a thrust fault ramp: numerical mechanical models. Journal of Structural Geology 23, 11–23. Ericsson, J.B., McKean H.C., Hooper R.J., 1998. Facies and curvature controlled 3D fracture models in Cretaceous carbonate reservoir, Arabian Gulf. In: Jones, G., Fisher, Q.J., Knipe, R.J. (Eds.), Faulting, Fault Sealing and Fluid Flow in Hydrocarbon Reservoirs, Special Publication 147. Geological Society, London, pp. 299–312. Jamison, W.R., 1997. Quantitative evaluation of fractures on monkshood anticline, a detachment fold in the foothills of Western Canada. American Association of Petroleum Geologists Bulletin 81, 1110–1132. Lisle, R.J., 1994. Detection of zones of abnormal strains in structures using Gaussian curvature analysis. American Association of Petroleum Geologists Bulletin 78, 1811–1819. Lisle, R.J., 2000. Predicting patterns of strain from threedimensional fold geometries: neutral surface folds and forced folds. In: Cosgrove J. M., Ameen M. S. (Eds.), Forced Folds and Fractures, Special Publication 169. Geological Society, London, pp. 213–221. Murray, G.H., 1968. Quantitative fracture studyFsanish pool mckenzie country. North Dakota. American Association of Petroleum Geologists Bulletin 52, 57–65. Salvini, F., Storti, F., 2001. The distribution of deformation in parallel fault related fold with migrating axial surfaces: comparison between fault-propagation and fault-bend folding. Journal of Structural Geology 23, 25–32. Schultz-Ela, D.D., Yeh, J., 1992. Predicting fracture permeability from bed curvature. In: Tillerson, J.R., Wawersik, W.R. (Eds.), Proceedings of the 33rd US Symposium on Rock Mechanics, pp. 579–589. Snow, D.T., 1969. Anisotropic permeability of fractured media. Water Resources Research 5 (6), 1273–1289. Stearns, D.W., 1968. Certain aspects of fracture in naturally deformed rocks. In: Reicher, R. E., (Ed.), National Science Foundation Advanced Science Seminar in Rock Mechanics. Air Force Cambridge Research Laboratory Special Report, pp. 97–118. Stewart, S.A., Podolski, R., 1998. Curvature analysis of gridded geological surfaces. In: Coward, M.P., Daltaban, T.S., Johnson, H. (Eds.), Structural Geology in Reservoir Characterization, Special Publication 127. Geological Society, London, pp. 133–147.
