Understanding gravity gradients—a tutorial

Understanding gravity gradients ― a tutorial AFIF H. SAAD, Saad GeoConsulting, Richmond, Texas, USA THE METER READER THE LEADING EDGE AUGUST 2006 THE METER READER Coordinated by Bob Van Nieuwenhuise Understanding gravity gradients—a tutorial AFIF H. SAAD, Saad GeoConsulting, Richmond, Texas, USA T he use of gravity gradient (GG) data in exploration is becoming more common. However, interpretation of gravity gradient data is not as easy as the familiar vertical gravity data. For a given source, regardless of its simplicity, gravity gradients often produce a complex pattern of anomalies (single, doublet, triplet, or quadruplet) as compared to the simple single (monopolar) gravity anomalies. This paper is a minitutorial on gravity gradients and is designed to provide a simple explanation of the complex pattern of GG anomalies and suggest some guidelines for the interpretation of measured surface GG data. To demonstrate the complex pattern of anomalies associated with gravity gradients, I will compute the gravity gradient components of the full gradient tensor starting with the basic building block, the gravitational potential. This will be followed by computing and examining: • the first derivatives of the potential in x, y, and z directions (i.e., the horizontal and vertical components of the gravity field vector) • the second derivatives of the potential (x-, y-, and zderivatives of each gravity vector component) which constitute the nine components of the full GG tensor (of which only five are independent). Figure 1 shows the model, constructed with GOCAD, used for these computations—a diapiric salt body in a sedimentary section whose density increases with depth, in a geologic setting typical of the U.S. Gulf coast. The upper part of the salt is above the “nil” zone and, thus, has positive density contrasts with the surrounding sediments; the lower part of the salt body has negative density contrasts. The nil zone, at depth of about 1 km in this example, is the area where the density of the surrounding sediments is identical to that of salt; hence, its gravity effect is nil (Figure 2). This model is very realistic and useful because it was digitized from a real case history and is really two models in one—a shallow one with positive density contrasts, and a deeper one with negative density contrasts. Hence, it is useful for testing the resolving capabilities of gravity gradients from shallow to deep sources. The gravitational potential and its first derivatives. Figure 3 shows color contour maps of the gravitational potential (P) and its first derivatives in the x, y, and z directions (P,x; P,y; and P,z). These derivatives are the horizontal (P,x; P,y) and vertical (P,z) gravity components of the gravity field vector. The salt model depth contours (0.2, 0.5, 1, 2, 3, 4, 5, and 6 km) are projected on all the maps for reference and to aid in interpretation. The potential (P) shows mainly a broad bell-shaped negative anomaly due to the main salt body; the effect of the shallow part of the salt is not obvious although, on closer examination, there is a subtle change in the contour spacing in the northeast, suggesting a small positive anomaly. It is interesting to note that, in spite of the apparent simplicity of the potential anomaly, it contains all the informa942 THE LEADING EDGE AUGUST 2006 Figure 1. GOCAD salt model and Cartesian coordinates system used. Figure 2. Density-depth curves for salt and sediments typical of Gulf of Mexico geologic setting. tion that produces the enhanced details and complex anomalies of the gravity and gravity gradient components shown later. The first horizontal derivatives of the potential in x and y or E and N directions produce doublet anomalies, a negative–positive pair along the x and y axes, respectively (Figure 3, top row). These are equivalent to the horizontal gravity components gx and gy that would be measured by Figure 3. Gravitational potential P and its first derivatives P,x, P,y, and P,z (x-, y-, and z-gravity field components of the gravity vector g due to the salt model shown). Figure 5. Frequency responses and characteristics of first derivative filters: horizontal derivatives (left), vertical derivative (right). Figure 4. First horizontal derivative of P in the NE direction. a horizontal gravimeter. The pattern of doublet anomalies is coordinate-dependent as suggested by the rotated pattern in Figure 4 for the NE directional horizontal derivative. We should expect this pattern of gravity anomalies if we consider the characteristic properties of the horizontal derivatives. The horizontal derivative operator is a phase filter (left panel in Figure 5) which will shift the location of anomalies or, in this case, split the negative P anomaly into a negativepositive pair along the x- or y-axis, respectively. The frequency response of Ꭿ/Ꭿx, for example, is ikx where i is the imaginary number, and kx is the wavenumber in the x direction. Hence, the x-derivative involves a phase transformaAUGUST 2006 THE LEADING EDGE 943 Figure 6. Gravity gradients (second derivatives of the potential). tion as well as enhancement of high frequencies (or high wavenumbers) relative to low frequencies. The phase transformation generally produces anomaly peaks (or troughs) approximately over the source edges in the case of wide bodies (width w is large relative to depth d, w > d). The enhancement of high wavenumbers sharpens these peaks to increase the definition of body edges in addition to emphasizing the effects of shallow sources. Another explanation, from elementary calculus, is that in the space domain the horizontal derivative is defined as the rate of change of P with respect to x or y. Hence, the horizontal derivative is a measure of the slope or “gradient” of the anomalies in the x or y direction (Figure 5, bottom left). If we consider the P surface as topography, 944 THE LEADING EDGE AUGUST 2006 Figure 7. Frequency responses of second derivative filters. the potential P (Figure 3) has a negative slope on the west and south sides of the minimum (going “downhill”), zero slope at the minimum, and positive slope on the east and north sides (going “uphill”)—thus producing the negativepositive pairs of gravity anomalies P,x and P,y. Notice that we can obtain the P,y pattern of anomalies by a simple 90° counterclockwise rotation of the P,x pattern, in the same manner as one rotates the x-axis to the y-axis. In fact, if we rotate the x- and y- axes 45° counterclockwise, or if we take the directional horizontal derivative of the potential P in the NE direction, the negative-positive pattern of anomalies obtained is rotated in the same direction as shown in Figure 4, emphasizing the fact that these anomalies are coordinate-dependent. Figure 8. The full gravity gradient tensor. The first vertical derivative, on the other hand, is a zerophase filter (right panel of Figure 5); hence, it will not affect the location of anomaly peaks, but it will sharpen the potential anomalies and will emphasize the high-frequency components due to shallow sources relative to the deeper effects, as seen in the P,z map of Figure 3 (lower right). The vertical derivative of P is, by definition, the rate of change of P with depth; hence, its effect will be similar to downward continuation, making the anomalies sharper and emphasizing shallower effects. Notice that the P,z data are the vertical gravity component gz measured by modern-day gravimeters. The frequency response of all three first derivative filters (Figure 5) is proportional to the wavenumber; hence, we expect these derivatives to enhance the short wavelengths or high frequencies due to the shallow part of the salt with positive density contrast as suggested by the bending or embayment of the contours at that location in Figure 3. Notice the reverse polarity of the shallow anomalies in response to the positive density contrast of the salt as compared to the deeper salt effect. The second derivatives of the potential. The various gravity gradient components are computed by taking the horizontal x- and y-derivatives and vertical z-derivative of each of the three gravity components of Figure 3. Figure 6 shows the five independent components of the gravity gradient tensor (second derivatives of the potential P): P,xx; P,xy; P,xz; P,yy; and P,yz along with the dependent second vertical derivative P,zz (P,zz = –P,xx –P,yy by Laplace’s equation). Again, we can expect the single, double, triple, and quadruple pattern of anomalies produced, if we keep in mind the properties and effect of the derivative operators explained above, or the frequency responses of the second derivative filters shown in Figure 7. The gravity gradient component P,xx is computed by taking the x-derivative of P,x. This results in a second phase AUGUST 2006 THE LEADING EDGE 945 Figure 9. Combined products of gravity gradient components: Horizontal gradient and total gradient of gz. Figure 10. Combined products of gravity gradient components: Differential curvature magnitude. 946 THE LEADING EDGE AUGUST 2006 Figure 11. Gravity gradient invariants (after Pedersen and Rasmussen, 1990). Figure 12. Other gravity gradient combinations: Euler deconvolution using GG tensor components. transformation and further enhancement of the high frequencies of the anomalies of P,x. Thus, the negative anomaly of the doublet of P,x splits into a negative-positive pair, and the positive anomaly splits into a positive-negative pair west-to-east along the x-axis, resulting into a “negativestrong positive–negative” triplet (P,xx of Figure 6). We can also explain this pattern by examining the slopes of the anomalies of P,x as we proceed from left to right along the x-axis. Notice that the steepest slope is at the center of the map of P,x (Figure 3) and it is positive; the zero slopes are at the trough and peak of P,x, and the gentle negative slopes are to the left and right of the trough and peak, respectively. In a similar manner, we can explain the triplet pattern of the component P,yy (center panel of Figure 6) which is simply a 90° counterclockwise rotation of the P,xx pattern. The gravity gradient component P,xy is computed by taking the derivative of P,x in the y (or N) direction or by taking the derivative of P,y in the x (or E) direction. This results in a second phase transformation and further enhancement of high-frequencies of the anomalies of P,x or P,y. Considering P,x, the negative anomaly of the doublet of P,x splits into a negative-positive pair, along the y-direction or south-to-north and the positive anomaly splits into a positive-negative pair along the y-direction or south-to-north, resulting in a “negative-positive–negative–positive” quadruplet (P,xy of Figure 6, top center panel). We can also explain this pattern by examining the slopes of the anomalies of P,x in Figure 3 as we proceed from south-to-north in the y-direction, or the slopes of the anomalies of P,y as we proceed from west-to-east in the x-direction. The gravity gradient components P,xz and P,yz and P,zz (right column of Figure 6) are computed by taking the z-derivative of P,x and P,y and P,z, respectively. This only causes further sharpening of the anomalies and enhancements of the high frequencies of P,x and P,y and P,z without any changes in the location or shapes of the anomalies, the z-derivative being a zero-phase filter (Figure 7). The full gradient tensor can be constructed by noting that P,yx = P,xy and P,zy = P,yz and P,zx = P,xz (Figure 8). The tensor is symmetric about its diagonal and its trace, the sum of the diagonal components (P,xx + P,yy + P,zz), is identically equal to zero in source-free regions, according to Laplace’s equa- tion. Thus, the tensor has only five independent components. It is interesting to note from Figure 8 that the first (top) row of the tensor is identical with the first (left) column and its components are the x-, y- and z-derivatives of the gravity field horizontal component gx of the gravity vector g (Figure 3). Similarly, the second (center) row of the tensor is identical with the second (center) column and its components are the x-, y- and z-derivatives of the horizontal gravity field component gy of the gravity vector g; the third (bottom) row of the tensor is identical with the third (right) column and its components are the x-, y- and z-derivatives of the gravity field vertical component gz of the gravity vector g. Notice the greater enhancement and better definition of the shallow anomaly pattern associated with the upper part of the salt in all gravity gradient maps (Figure 6). This is because the frequency response of all second derivative filters is proportional to the square of the wave number (Figure 7). Notice also the reverse polarity of the high-frequency anomaly pattern in all components as expected from the positive density contrast of the shallow salt. Thus, for example, the triplet of P,xx is “positive-negative-positive” for shallow salt as compared to the main “negative-positive-negative” pattern for the deep salt. One should emphasize that the pattern of anomalies produced is coordinate-dependent. However, one can use these patterns and shapes of gravity gradient anomalies with the projected outline of the causative salt body in this example to develop interpretation techniques for locating the main salt body, its edges, and its shallow part. For example, the zero contours of P,xx and P,yy closely define the westeast edges and south-north edges of the main salt body, AUGUST 2006 THE LEADING EDGE 947 tions discussed above do not hold in this case; the locations of the zero contours, the lows and highs, and size of the anomalies in general will depend mainly on the depth to the source, rather than the width/depth ratio. Finally, the P,zz anomalies (Figure 6, bottom-right panel) can be used to locate the center of the anomalous source mass. Combinations of GG components (invariants). Various combinations of the gravity gradient components can be used to simplify their complex pattern and to further enhance and aid in the interpretation of the data. Figures 9 and 10 show three examples: amplitude of the horizontal gradient of vertical gravity (gz); amplitude of the total gradient or analytic signal of gz; and the differential curvature which is also known from the early torsion balance literature as the horizontal directive Figure 13. Similarity between surface horizontal gravity (in the X-Y plane) and subsurface vertical tendency or HDT. The horizontal and gravity (in the X-Z plane). total gradients of gz (Figure 9) are computed from combinations of the elements of the third column (or third row) of the gravity gradient tensor— P,xz and P,yz and P,zz (Figure 6). The latter are the x, y, and z derivatives of P,z (or gz). The horizontal gradient of gz can be used as an edge-detector or to map body outlines. The analytic signal can be used for depth interpretation. The differential curvature (Figure 10) is computed by a combination of the other components of the tensor: P,xx and P,xy and P,yy. The magnitude of the differential curvature emphasizes greatly the effects of the shallower sources. Several interpretation techniques for the differential curvature are available in the early literature of the torsion balance. The three examples of combined GG products discussed above are useful in simplifying and “focusing” the complex pattern of anomalies over their source, providing more enhancements to the high-frequency part of Figure 14. Similarity between surface horizontal gravity gradient difference (in the X-Y plane) and anomalies due to shallow sources, and subsurface vertical gravity gradient (in the X-Z plane). producing coordinate-independent or invariant anomalies. These are perrespectively (Figure 6, top left and center panels). Also, the haps easier to interpret than the original gradient compopeaks and troughs of the quadruplet pattern of P,xy anom- nents. Other coordinates-independent invariants can be alies are located roughly around the perimeter of the salt computed and used as well for interpreting the data using body (Figure 6, top center) and can be used to delineate the different combinations of the GG components. For examsalt boundary. The negative-positive pairs of the P,xz and ple, one can compute the horizontal and total gradients of P,yz anomalies are near or on the west-east and south–north gx and gy from the elements of the first row and second row edges of the body, respectively (Figure 6, top-right and cen- of the GG tensor, respectively. Figure 11 defines other gravter-right panels). These relations depend on the width/depth ity gradient invariants, I0, I1, and I2 suggested by Pedersen ratio of the source and are generally valid only for wide bod- and Rasmussen (1990) and used for interpretation of GG ies, i.e., bodies whose width is greater than their depth (w> data. Gravity gradient components can also be combined d). It should be emphasized that narrow sources (w≤d), to form three different Euler equations for gx, gy, and gz that including point masses, will produce similar geometric pat- can be used to solve for source depth (Figure 12), as sugtern of complex anomalies as in Figure 6; however, the rela- gested by Zhang et al. (2000). 948 THE LEADING EDGE AUGUST 2006 tance X=2R from the center of the sphere of radius R. The apparent gravity doublet and GG triplet patterns encountered in the borehole are similar to the patterns of gravity gradient profiles that would be observed on the surface. Thus, interpretation techniques developed and used for borehole gravity and gravity gradient data can be extended and used for surface gravity gradient data interpretation. Overall, experience with interpretation of borehole gravity data can be valuable for the interpretation of surface gravity gradient profile and map data. Figure 15. Borehole vertical gravity and gravity gradient (apparent density) profiles due to a sphere of radius R, density contrast ∆ρ. Borehole distance X = 2R from the center of the sphere. Similarities between surface and subsurface gravity and gravity gradients. It is interesting to note that there are similarities between surface variations of the horizontal gravity and GG components and subsurface variations of vertical gravity and vertical GG (or anomalous apparent density) such as those observed in a borehole. Figures 13 and 14 show examples illustrating these similarities. Figure 13 compares surface variations in the x-y plane of the horizontal gravity component P,y (Figure 3) with subsurface variations in the x-z plane of vertical gravity due to a spherical source. Figure 14 shows a similar comparison between surface gravity gradient difference (P,xx – P,yy) and subsurface vertical gravity gradient or apparent density anomaly, as used in borehole gravity work, due to the same spherical mass. Vertical profiles in the z direction extracted from the maps on the right-hand sides of Figures 13 and 14 show the anomalous responses expected in boreholes and measured in borehole gravity surveys (Figure 15). In this example, the boreholes are located at a remote dis- Conclusions. Gravity gradients (GG) often produce a pattern of complex anomalies that is coordinate-dependent, not necessarily reflecting the shape of the underlying sources. Understanding GG anomalies is important in the interpretation of measured data. It is easy to understand the complex pattern of gravity gradients if one considers the fact that they are derivable from the simple gravitational potential, being the directional second derivatives of the potential. In general, for 3D sources producing single bell-shaped potential and vertical gravity anomalies, the P,zz gravity gradient component consists of a single anomaly; the P,xz and P,yz components consist of doublet anomalies; the P,xx and P,yy components consist of triplet anomalies; and the P,xy component consists of quadruplet anomalies. Various combinations of GG components can be used to produce coordinate-independent “invariants” that are simple, easy to interpret, more localized, and more related to the size and shape of the sources. There are also similarities between surface and subsurface (or borehole) variations of certain gravity and gravity gradient components. Hence, interpretation methods developed and used for borehole gravity data may be applicable or can be extended to surface GG data interpretation. Certainly past experience with borehole gravity can be valuable in interpreting surface gravity gradient data. Suggested reading. “Gravity gradiometry resurfaces” by Bell et al. (TLE, 1997). “Gravity gradiometry in resource exploration” by Pawlowski (TLE, 1998). “The gradient tensor of potential field anomalies: Some implications on data collection and data processing of maps” by Pedersen and Rasmussen (GEOPHYSICS, 1990). “Euler deconvolution of gravity tensor gradient data” by Zhang et al. (GEOPHYSICS, 2000). TLE Acknowledgments: Parts of this work were conducted while the author was employed by Gulf Research and Development, Chevron, and Unocal companies. This paper was presented at the SEG75 Annual Meeting in Houston, Texas. Corresponding author: afifhsaad@netscape.net AUGUST 2006 THE LEADING EDGE 949 Author Biography Afif H. Saad is a Geophysical Consultant, specializing in integrated gravity / magnetic / seismic / geologic interpretation, modeling, magnetic depth estimation, software development, and training. He has over 25 years of experience in the oil industry, including GULF R&D, GULF E&P, CHEVRON and UNOCAL Oil Companies. He also held positions with Aero Service Corp. in Philadelphia and LCT Inc. in Houston as well as in the academia at Cairo University, Stanford University, and University of Missouri at Rolla. Afif received a Ph.D. in Geophysics from Stanford University, M.S. in Geology/Geophysics from Missouri School of Mines-Rolla, and B.Sc. (Honors) Special Geology from Alexandria University, Egypt. He is a member of SEG, Gravity and Magnetics Committee, and GSH. He was the chairman of the Houston Potential Fields SIG of GSH from 2000-2004, and an Associate Editor for GEOPHYSICS – Magnetic Exploration Methods from 1999 – 2005.