Application of the finite element method to gravity data case study: Western Turkey

Journal of Geodynamics 39 (2005) 431–443

Application of the finite element method to gravity data case study: Western Turkey Ilknur Kaftan ∗ , Mujgan Salk, Coskun Sari Dokuz Eylul University, Faculty of Engineering, Departmant of Geophysics, 35160 Buca-Izmir, Turkey Accepted 18 April 2005

Abstract Gravity anomalies always include the total effects (combination of the structures which have different densities and depths) of the study area and beyond. And the well-known non-uniqueness of potential field modelling may lead to very different interpretation results. The finite element method (FEM), which has been used in potential field interpretation for decades, makes complex problems to be solved easily and accurately. The first step of FEM is to identify the elements and then to decide on the boundary of the solution space. In this step, the solution space is divided into elements. After determination of the geometrical structure of the solution space, the most suitable elements should be chosen for this geometrical structure. The agreement between the geometry and the elements is quite important for the convergence to the best possible solution. In this work, the methods of trend analysis, filtering, analytical continuation and FEM were applied to a theoretical model and to gravity data from western Turkey to produce the regional and the residual anomalies. The results were compared and it was found that the FEM produced more accurate results than other methods did. © 2005 Elsevier Ltd. All rights reserved. Keywords: Finite element method; Regional/residual separation; Comparison of methods

1. Introduction Geophysical data processing methods, such as trend analysis, analytical continuation and filtering can reveal the general structural properties of a region. Although, there were improvements in distinguishing ∗

Corresponding author. Tel.: +90 232 4531008/1207; fax: +90 232 4538366. E-mail address: [email protected] (I. Kaftan).

0264-3707/$ – see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jog.2005.04.003

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

432

the regional and the residual gravity anomaly, no progress has been made in trend analysis, filtering or analytical continuation methods. By the improvements of scientific computing, high degree trend analysis and filtering methods have been used more efficiently for 2-D data. On the other hand, computed regional anomalies still contain residual anomaly effects. For this reason, finite element method (FEM), might be used to eliminate the residual anomaly effects on the regional anomalies. As much as the potential theory may allow. The gravity map in the real space is superimposed to compute the regional anomalies on the finite elements in the FEM application. For each node coordinate and field variable, the element can be defined by the same shape functions. Due to this feature, it could be named an isoparametric element. The observed gravity values for each node are obtained from the map. The non-dimensional reference plane, which is related to the real plane, was described by using the shape function to perform an easy computation. Eight observed gravity values on quadratic isoparametric element, which superimposed the gravity map are needed to compute regional. Other observed gravity data are not needed from the gravity map. The results from FEM, trend analysis, filtering and analytical continuation methods are compared and it is seen that FEM was capable of separating the regional and residual anomaly more accurately than the traditional methods were.

2. The computation of gravity data Observed gravity data are the total effect of deep and shallow structures (Pawlowski, 1994) and can be written as g(x, y) = gs (x, y) + gd (x, y)

(1)

where, g(x,y) are the observed gravity data, gs (x,y) the gravity effect of shallow structures and gd (x,y) is the gravity effect of deep structures. One of the problems in the gravity interpretation is to distinguish the regional and the residual anomalies. It may provide a more accurate definition of the structures. There are several techniques to compute the regional anomalies from the observed gravity data, such as trend analysis, filtering and analytical continuation. Although there were some improvements in these techniques in the past decades, it could not be possible to separate residual effects from the regional field exactly due the intrinsic property of potential theory. As mentioned earlier, FEM is here applied in order to determine the regional anomalies. The method has advantages over the traditional methods. For instance, only a few observation points on a Bouguer gravity map are needed to compute the regional gravity anomaly. The FEM can be applied for any size of gravity map. In addition, it will be shown that the computed regional anomaly contains minimal residual anomaly effects.

3. Theory of FEM FEM has been successfully applied to different modelling problems in the scientific world including in geophysics for decades. In this study, it was used to distinguish the regional and residual anomalies by using the necessary shape functions. These shape functions play quite an important role to change the

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

433

Fig. 1. Eight-node quadrilateral element in the x–y plane (taken from Mallick and Sharma, 1999).

reference plane. After the method was applied, it was seen that it produced better results than the existing traditional methods produced. In addition to that, fewer gravity values were needed. For instance, only eight nodes with the shape functions were used to compute the regional anomaly. The regional anomaly was computed (Mallick and Sarma, 1992; Sarma et al., 1993) and given as g(x, y) =

8


Ni (x, y)gi

(2)

i=1

where g(x,y) represents the regional field, Ni (x,y) represents the shape functions and gi represent the node values. An eight node quadrilateral element in the (x,y) plane is shown in Fig. 1, where the real map space (xc ,yc ) represents the center of the eight node quadrilateral element. The length of the sides are 2a and 2b. An eight node quadrilateral element in the non-dimensional (ξ,η) plane is shown in Fig. 2. This new reference plane was used to compute the shape functions. As is illustrated in Fig. 2, changing the existing (x,y) plane coordinates to (ξ,η) plane coordinates is a necessary step.

Fig. 2. Eight-node quadrilateral element in the non-dimensional ξ,η plane (taken from Mallick and Sharma, 1999).

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

434

Fig. 3. Triple cubic blocks model.

The new (ξ,η) plane coordinates can be expressed (Mallick and Sharma, 1999) as x − xc y − yc η= a b The shape functions from 1 to 8 nodes in Fig. 2 can be expressed (Cheung and Yeo, 1979) as ξ=

Ni (ξ, η) =

(1 + ξξi )(1 + ηηi )(ξξi + ηηi − 1) 4

Fig. 4. Theoretical gravity anomaly map (contour interval is 0.5 mGal).

(3)

(4a)

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

435

for i = 1, 3, 5 and 7 only (1 − ξ 2 )(1 + ηηi ) 2 for i = 2 and 6 only Ni (ξ, η) =

(1 + ξξi )(1 − η2 ) 2 for i = 4 and 8 only It should be noted that Ni (ξ,η) = 1 at the ith node and zero at the other nodes, therefore Ni (ξ, η) =




(4b)

(4c)

Ni (ξ, η) = 1

The regional gravity could be expressed in terms of the shape functions (Mallick and Sharma, 1999) as g(ξ, η) =

8


Ni (ξ, η)gi

(5)

i=1

Fig. 5. x–y plane on the theoretical gravity anomaly map.

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

436

where gi represents the eight nodes gravity values on the gravity map and Ni (ξ,η) represents the shape functions. Based on the new reference plane, the new coordinates could be given (Mallick and Sharma, 1999) as x(ξ, η) =

8


Mi (ξ, η)xi

(6a)

i=1

Fig. 6. (a) Upward continuation anomaly map for h = 10 km (contour interval is 0.5 mGal). (b) Low-pass filtering anomaly map (cut-off frequency is 0.1 cycle/grid spacing and contour interval is 0.5 mGal). (c) The third order trend analysis map (contour interval is 0.5 mGal). (d) Finite element regional anomaly map (contour interval is 0.5 mGal).

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

y(ξ, η) =

8


Mi (ξ, η)yi

437

(6b)

i=1

where Mi (ξ,η) represents the shape functions, xi and yi represent the new node coordinates (Fig. 3).

Fig. 7. (a) Downward continuation anomaly map for h = 1 km (contour interval is 0.5 mGal). (b) High-pass filtering map (cut-off frequency is 0.1 cycle/grid spacing and contour interval is 0.2 mGal). (c) Finite element residual anomaly map (contour interval is 0.5 mGal).

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

438

These equations provide the transformation of the computed regional anomaly from the reference (ξ,η) plane to the real (x,y) plane. In Eqs. (5), (6a) and (6b), it can be seen clearly that the coordinates of any given point and field variable use the same shape functions, that is, Ni (ξ,η) = Mi (ξ,η), called isoparametric elements. This is one of the important features of the isoparametric elements. 4. Application of FEM FEM was applied to both the theoretical gravity anomaly and to a Bouguer gravity anomaly from Western Turkey (Kaftan, 2003). As shown in Fig. 5, an eight nodes quadrilateral element was superimposed on the theoretical gravity map. For each node, a gravity value was obtained from the gravity map. Each gravity value at these nodes, say 1–8, was transferred to the reference (ξ,η) plane. Then, the regional anomaly was computed for each node by using Eq. (2) in the reference (ξ,η) plane. In this step, only the values of the shape functions got changed for each node in the (ξ,η) reference plane. There was no change in the gravity values obtained from the gravity anomaly map. Once this step was completed, it was possible to change the coordinates back to the (x,y) plane by using Eqs. (6a) and (6b) and then to compute the regional anomaly. Finally, the residual anomaly was computed by using Eq. (1). 4.1. Theoretical model Three cubic blocks with different volumes as shown in Fig. 3 were used to compute the theoretical gravity anomaly and the computed anomaly. The 3-D gravity anomaly of any given location is expressed

Fig. 8. The Bouguer gravity anomaly map from Western Turkey (taken from MTA, 1979, contour interval is 10 mGal).

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

439

Fig. 9. (a) The fifth order trend analysis anomaly map of Western Turkey (contour interval is 10 mGal). (b) Upward continuation anomaly map from Western Turkey for h = 1 km (contour interval is 7 mGal). (c) Low-pass filtered anomaly map from Western Turkey (cut-off frequency is 0.1 cycle/grid spacing and contour interval is 10 mGal). (d) Finite element method regional anomaly map from Western Turkey (contour interval is 10 mGal).

440

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

Fig. 10. (a) Downward continuation anomaly map from Western Turkey for h = 1 km (contour interval is 60 mGal). (b) High-pass filtering anomaly map from Western Turkey (cut-off frequency is 0.1 cycle/grid spacing and contour interval is 5 mGal). (c) Finite element method residual anomaly map from Western Turkey (contour interval is 5 mGal).

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

441

analytically (Grant and West, 1965) as, B00 z0 + 3B20 [D + E − F + G − H] /2R + 3B22 [I + K + L − M + N] /R7 (7) R3 where, a, b, c,ρ, z, w and Ψ represent the half width in the x-direction, half length in the y-direction, half thickness in the z-direction, density contrast, depth of the cubic block, the angle between the cubic block and the x-axis, the angle between the cubic block and the z-axis, respectively. Using Eq. (7) (see Appendix A), the computed theoretical gravity anomaly map is shown in Figs. 4 and 5. The regional anomaly map for the theoretical gravity anomalies is shown in Fig. 6. As shown in Fig. 6(d), the regional effect of the third deep cubic block are more visible in the FEM application. On the other hand, trend analysis, downward continuation and low-pass filtering were applied to the theoretical anomaly map with are shown in Fig. 7. The shallow structure effects can be seen clearly on both the FEM application and also on the traditional method applications. g(x, y) =

4.2. Bouguer gravity map of Western Turkey The FEM and the other traditional methods were applied to the Bouguer gravity data from Western Turkey to compute the regional and the residual anomalies. The study area is between 38◦ and 40◦ latitude and between 28◦ and 29◦ longitude, namely in the Aegean graben systems (Fig. 8). The regional anomaly maps from the Bouguer gravity map of Western Turkey are shown in Fig. 9. From Fig. 9(a–c), it is clearly seen that the residual anomaly affects the regional anomaly in the trend analysis, in upward continuation and low-pass filtering applications, but, as shown in Fig. 9(d), the FEM regional anomaly map contains less residual anomaly effects than the other traditional methods have. The residual anomaly maps from the Bouguer gravity map of Western Turkey are shown in Fig. 10. Due to the structural effects of the target, the shallow structures can not be seen in the upward continuation application results, as illustrated in Fig. 10(a), but the Aegean graben systems can be seen on the high-pass filtering and the FEM application results (Fig. 10b and c).

5. Conclusion Because of the non-uniqueness nature of the problem in the solution space, gravity method has always some problems in the interpretation stage. Trend analysis, filtering and analytical continuation methods have been applied to eliminate these problems for many years. In this study, the trend analysis, filtering, analytical continuation and the FEM were applied to both the Bouguer gravity data from Western Turkey and the theoretical gravity data. The finite element method superimposed the gravity anomaly map to compute the regional anomaly by using the shape functions. When all the results were compared, it was seen that the finite element method was capable of separating the regional and the residual anomalies more accurately and effectively than the traditional methods were. Especially, in the synthetic data application, it was seen that the regional anomaly had less residual anomaly effects on it. It can be clearly seen that the problem of regional–residual separation is one of interpretation and thus subject to the principal and unsolvable ambiguity of gravity inversion. No method, however sophisticated, can give us a regional or a residual field, which are unambiguously due to unique mass distributions.

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

442

Nevertheless, different methods perform differently and their geophysical–geological value depends on the geological problem at hand.

Acknowledgements We would like to thank Dr. K. Mallick and K.K. Sharma about applying the finite element methods and ¨ for helping with the computer code. We also would like to thank Dr. Zafer Y¨ucel Oner and Prof. Dr. Wolf Jacoby for reviewing our manuscript. This research project was funded by the Dokuz Eylul University Research Fund, Project Code: 03.KB.FEN.026.

Appendix A The terms of Eq. (7) can be expressed below, B00 = 8Gρabc; B20 =

B00 (2c2 − a2 − b2 ) 6

B22 =

B00 (a2 − b2 ) 12

D = (3 cos2 ψ − 1)z30 + (5 sin2 w sin2 ψ − 2 cos2 ψ − 1)x2 z0 E = (5 cos2 w sin2 ψ − 2 cos2 ψ − 1)y2 z0 F = 2 sinw sinψ cosψ(x3 + xy2 − 4xz20 ) G = 2 cosw sinψ cosψ(y3 + x2 y − 4yz20 ) H = 10 sinw cosψ sin2 ψ(xyz0 ) I = −3 sin2 ψz30 + (5 cos2 w − 5 sin2 w cos2 ψ + 2 sin2 ψ)x2 z0 K = (5 sin2 w − 5 cos2 w cos2 ψ + 2 sin2 ψ)y2 z0 L = 2 cosw sinψ cosψ(y3 + x2 y − 4yz20 ) M = 2 sinw sinψ cosψ(x3 + xy2 − 4xz20 ) N = 10 sinw cosw(1 + cos2 ψ)xyz0 1/2

R = (x2 + z20 )

I. Kaftan et al. / Journal of Geodynamics 39 (2005) 431–443

443

References Cheung, Y.K., Yeo, M.F., 1979. A Practical Introduction to Finite Element Analysis. Pitman Publishing, Marshfield, MA, p. 180. Grant, F.S., West G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill Book Company, New York, St. Louis, San Francisco, Toronto, London, Sydney. Kaftan, I., 2003. Application of Finite Element Method on gravity data in Western Anatolia, Dokuz Eylul University, The Graduate School of Natural and Applied Sciences (Master Thesis, in Turkish). Maden Tetkik Arama, M.T.A, Directorate of Mineral Research and Exploration, Turkey, 1979. The Bouguer Gravity Map from Western Turkey, unpublished. Mallick, K., Sarma, G.S., 1992. Finite element formulation for seismicwave propagation in oil-bearing geological formations. In: IMACS International Symposium on Mathematical Modelling and Scientific Computation, pp. 1–12. Mallick, K., Sharma, K.K., 1999. A finite element method for computation of the regional gravity anomaly. Geophysics 64 (2), 461–469. Pawlowski, R.S., 1994. Green’s equivalent-layer concept in gravity band-pass filter design. Geophysics 59, 69–76. Sarma, G.S., Gadhinglajkar, V.R., Mallick, K., 1993. Finite element simulation of bright-spot structures. J. Assoc. Expl. Geophys. 14, 43–47.