3D inversion of gravity data by separation of sources and the method of local corrections: Kolarovo gravity high case study

Journal of Applied Geophysics 75 (2011) 472–478

Contents lists available at SciVerse ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

3D inversion of gravity data by separation of sources and the method of local corrections: Kolarovo gravity high case study Ilya Prutkin a,⁎, Peter Vajda b, Robert Tenzer c, Miroslav Bielik b, d a

Institute of Earth Sciences, Jena University, Burgweg 5 11, 07749 Jena, Germany Geophysical Institute, Slovak Academy of Sciences, Dubravska cesta 9, 845 28 Bratislava, Slovakia National School of Surveying, Division of Sciences, University of Otago, 310 Castle Street, Dunedin, New Zealand d Department of Applied and Environmental Geophysics, Faculty of Natural Sciences, Comenius University, Mlynska dolina, 842 15 Bratislava, Slovakia b c

a r t i c l e

i n f o

Article history: Received 30 March 2011 Accepted 12 August 2011 Available online 14 September 2011 Keywords: Potential field methods 3D inversion Gravity interpretation Anomalous body

a b s t r a c t We present a novel methodology for 3D gravity/magnetic data inversion. It combines two algorithms for preliminary separation of sources and an original approach to 3D inverse problem solution. The first algorithm is designed to separate sources in depth and to remove the shallow ones. It is based on subsequent upward and downward data continuation. For separation in the lateral sense, we approximate the given observed data by the field of several 3D line segments. For potential field data inversion we apply a new method of local corrections. The method is efficient and does not require trial-and-error forward modeling. It allows retrieving unknown 3D geometry of anomalous objects in terms of restricted bodies of arbitrary shape and contact surfaces. For restricted objects, we apply new integral equations of gravity and magnetic inverse problems. All steps of our methodology are demonstrated on the Kolarovo gravity anomaly in the Danube Basin of Slovakia. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Our main goal was to develop an approach for interpretation of potential field data, which has to face the following challenges: i) to be really three-dimensional (to avoid an assumption that the object sought is an infinite cylinder); ii) to abandon the simplest bodies like prisms, ellipsoids, cylinders of finite length for objects of more complex geometry; iii) to take into account the non-uniqueness of the inverse problem. Our approach provides different admissible solutions to the geological section based on potential field data. Later it can be discriminated among these solutions based on additional geophysical, geological or tectonic constraints (if such information exists). We clarify all details of the approach while interpreting the Kolarovo gravity high data in the southern Slovak part of the Pannonian Basin. This gravity anomaly is located in the south-eastern part of the Danube Basin (a part of the Pannonian Basin in Central Europe). The anomaly belongs to one of the largest and most famous gravity highs in the Western Carpathian area. Therefore it has been of significant interest to geophysicists and geologists since early sixties of the twentieth century. In terms of the Bouguer anomaly, the Kolarovo high reaches a magnitude of 28 mGal. It is supposedly associated with deep seated sources within the pre-Neogene basement (Kubeš et al., 2010). According to Kubeš et al. (2010), the interpretations of

⁎ Corresponding author. Tel.: + 49 3641 948751; fax: + 49 3641 948662. E-mail address: [email protected] (I. Prutkin). 0926-9851/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2011.08.012

these sources are insufficient and problematic and thus remain open for further investigation. The quantitative interpretation of this gravity high (Bielik et al., 1986; Sitárová et al., 1984) has revealed that the density of the anomalous body varies from 2900 to 3050 kg/m3, compared to the density of the upper crust of 2700 kg/m3. Sitárová et al. (1984) determined a source body as a right-rectangular homogeneous prism with a top plane at 4.5 km with a center of mass at 9.5 to 12.5 km. Bielik et al. (1986) interpreted the Kolarovo anomaly by means of 3D gravimetric forward modeling using polyhedrons. The TFM (Truncation Filtering Methodology) interpretation of the Kolarovo high yielded an estimate of the depth of the center of mass of the body at 8.7 km (Vajda et al., 2002). We start with separation of sources. First, we remove shallow sources down to the depth corresponding to the assumed basement. Our algorithm is based on the subsequent upward and downward continuation of the observed gravity data. During this procedure, we integrate data along the area of investigation only. This has been made possible by prior subtraction of a model of the regional field, which we regard as a 2D harmonic function, from the observed data, thus forcing the data outside the investigation area to be zero. To illustrate the regional field subtraction and the separation of sources in depth, we refer to (Prutkin and Casten, 2009; Prutkin and Saleh, 2009), where this algorithm was applied to isolate the contribution of the Moho interface to the total field. Then, the 3D Moho topography was found by means of the method of local corrections applied to the thus separated residual field (Prutkin and Saleh, 2009). After obtaining the separated signal of interest, it is interpreted by means of direct gravity inversion, using the method of local corrections.

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

We do not use the time-consuming forward modeling approach. In the method of local corrections, a solution is calculated from gravity data automatically by successive approximations. The non-uniqueness and instability of the inverse problem are taken into account. Now we expand our methodology to be able to retrieve unknown geometry both for a contact surface and for a 3D restricted body of arbitrary shape. We approximate the given data by the field of several 3D line segments. This procedure allows to separate sources both in depth and also in the lateral direction, and to estimate their masses. In the next step it is possible to find a shape of a restricted object by means of the method of local corrections using new integral equations of gravity (and magnetic) inverse problems that were described in Prutkin (2008) and applied to the core–mantle boundary recovery. The approximation of the Kolarovo gravity anomaly with 3D line segment sources permits to study different models of the geological section. We are able to investigate various options to join segments and to relate the obtained gravitational effect to the same object (a restricted body or a contact surface). We have obtained a solution to the geological section which is in good agreement with results of previous interpretations mentioned above. Moreover, our new interpretation methodology provides three new admissible models: a body with more complicated geometry, an uplift of a relatively denser magmatic material as a contact surface, and a combination of an anomalous body and an elevation of the contact surface. 2. Separation of sources in depth We begin with the algorithm based on the upward and subsequent downward continuation (Prutkin and Casten, 2009), see also (Martyshko and Koksharov, 2005; Martyshko and Prutkin, 2003). The main purpose of this algorithm is to find a part of the observed field which is harmonic above a given depth h. We can treat this function as an effect of the half-space below the depth h. To find such function means to eliminate all sources located in the layer from the Earth’s surface, which we approximate by the plane z = 0, down to the prescribed depth. The algorithm allows separating the effects of shallow and deeper source objects. Firstly, we continue the data upward to the level z = h to diminish the influence of the sources in the near-surface layer. This operation causes the main errors in the vicinity of the boundary of the area. To reduce these errors we need a model of the regional field to be subtracted from the observed field prior to the upward continuation. Secondly, we continue the obtained function downward to the plane z = −h, i.e. over the distance 2h downwards. Since the downward continuation procedure is a linear ill-posed inverse problem, the regularization is applied. The function, which we treat as a regional field, is assumed to satisfy two-dimensional Laplace equation in the area of investigation S and to have the same values at the boundary of the area, as the observed field: 8 2 2 <∂ f ∂ f þ 2 ¼0 within area S 2 :∂ x ∂ y f ¼ Δg on its boundary ∂S

ð1Þ

If we subtract the values of this function, the residual field will be equal to zero at the boundary of the area. In this way no errors are introduced, when we integrate the residual field while upward continuation within the restricted area. According to the properties of harmonic functions, this function has no extreme values within the area, so we do not create artificial anomalies. Besides, as it is known from calculus of variations, a solution to the problem (1) minimizes the norm of the gradient, therefore, we obtain the smoothest possible function with given values on the boundary. All these properties allow us to regard this function as a model of the regional field.

473

After subtracting the suggested model of the regional field we perform upward continuation of the residual field by means of the Poisson integral:   1 h ∫∫  3=2 U ðx; y; 0Þdxdy ¼ U x′ ; y′ ; h    2 2 2π x−x′ þ y−y′ þ h2

ð2Þ

Formula (2) gives a solution of the Dirichlet problem: to find a harmonic function in the upper half-space with given values on the boundary (on the plane z = 0 ). At the second step we continue the obtained function downward to the depth h , i.e., over the distance 2h downwards. To do this, we apply a formula similar to Eq. (2):   1 h1 ∫∫  3=2 U ðx; y; −hÞdxdy ¼ U x′ ; y′ ; h ; ð3Þ    2π ′ 2 ′ 2 2 x−x þ y−y þ h 1

where h1 = 2h. Eq. (3) provides a function, which is harmonic in a half-space above the plane z = − h . This time we have to, and we do, treat the formula as an integral equation: the right hand side is given, while the unknown function U(x, y,−h) on the left hand side is sought. Although the signal of the shallow sources is attenuated after applying the upward continuation, it is not completely removed. As a consequence, the downward continuation is performed through the sources. It is worth emphasizing that the problem of downward continuation is a linear ill-posed inverse problem, therefore some regularization must be applied. Since the integral operator A in Eq. (3) is self-adjoint and positive, we apply the Lavrent'ev's approach (Lavrent'ev et al., 1986). If we write Eq. (3) in the following form: Au= uh, u is an unknown field at the level z = −h, uh is the obtained field at height h, then the Lavrent’ev’s regularization gives: (A + αI)u = uh, where I is the identity operator, and α is a regularization parameter. The choice of the regularization parameter decides how well the solution fits the data. To find a suitable regularization parameter, we use the method of residuals (Lavrent'ev et al., 1986), which exploits information about data noise: we look for the smoothest possible solution, which approximates the given data with the same accuracy as the data noise level. At last, we calculate (by an upward continuation) the field on the Earth’s surface z = 0 using a formula, similar to Eq. (3). We obtain a part of the field, which is harmonic down to the depth h, so that we now can treat it as an effect of the deeper sources, deeper than z =−h. We apply the described methodology to the Kolarovo gravity anomaly. Our goal is to remove shallow sources, thus we take for height h a value equal to 2 km, which corresponds to the minimum thickness of the Neogene sediments. The initial data and the results of our processing are shown in Fig. 1. The amplitude has practically not changed. The obtained field is very similar to the original one. It is a confirmation that the Kolarovo gravity high is caused mainly by deeper sources. 3. Approximation with 3D line segments To separate sources not in depth only, but also in the lateral direction, as well as to estimate their masses, we propose to approximate data by simple sources prior to 3D inversion. As a source for gravity data approximation, we use 3D line segments. The effect of a line segment can be evaluated by a quite simple formula. Let r0 = (x0, y0, z0) be the radius-vector of an observation point, r1 = (x1, y1, z1) and r2 = (x2, y2, z2)—radius-vectors of endpoints of a 3D line segment, w—its line density. We introduce following notations: α is the angle between vectors r0−r1 and r2−r1, c = |r2−r1| is the length of the line segment, c0 = |r0−r1|, c2 = |r0−r2|, a0 = z0−z1, a1 =

474

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

20 0

5

10

km

10

0

10 0

15

5

−5

−10

−20 −20

−30

−10

0

10

20

km −15 −11 −7

−3

1

9

5

13

17

21

25

mGal 20

0 5

10

km

10

0

0

10

5

15

−5

−10

−20 −30

−20

−10

0

10

20

km −15 −11 −7

−3

1

5

9

13

17

21

25

mGal Fig. 1. Initial data and their preliminary processing. Top-gridded data of the Kolarovo gravity high; bottom—contribution of shallow sources down to the basement is removed by means of the algorithm based on upward and downward continuation.

(z2−z1)/c. For the vertical derivative of the gravitational potential we have "  Vz ðr0 Þ ¼ w⋅ a1

1 1 − c0 c2

 þ

 # a1 c0 cos α−a0 c−c0 cos α   þ cos α ð4Þ c2 c20 1−cos2 α

for the general case cos2 α ≠ 1 (the observation point does not belong to the continuation of the line segment). In the special case cosα = 1 instead of Eq. (4) we obtain Vz ðr0 Þ ¼

" # w a0 þ a1 c0 a0 þ a1 ðc0 −2cÞ − 2⋅ c22 c20

ð5Þ

For cos α =−1 we have the same as in Eq. (5), but the endpoints of the segment change their places: c0 = |r0− r2|, c2 = |r0− r1|, a0 = z0−z2, a1 = (z1− z2)/c. The computation time to calculate the effect of a 3D line segment in one observation point exceeds the corresponding time for a point source only about 2 times. At the same time, it is a considerably flexible tool for the approximation of measurements: one 3D line segment approximates perfectly an effect of an object stretched in arbitrary direction; a short

segment generates a field similar to a point source, two segments are usually sufficient to approximate an effect of a restricted object with complex geometry. Besides, each 3D line segment is described by 7 parameters only: 3 coordinates of both two ends and the line density. This is an important point, because we obtain these parameters while approximating the data by means of nonlinear minimization. To demonstrate the advantages of the 3D line segment approximation, we begin with a synthetic field example. The gravitational effect of two rectangular prisms is calculated, one of them being cubic, the other being stretched in the vertical direction. In Fig. 2 we present the synthetic field (Fig. 2a) and projections of the prisms, as well as the obtained 3D line segments to the horizontal (Fig. 2b) and vertical (Fig. 2c) coordinate planes. Differences between the given data and the field of line segments (Fig. 2d) have RMS of 0.0188 mGal. This approximation yields not only the separation of sources and reliable mass estimates, but also reasonable information about the shape of the sought source body. Then we approximate the Kolarovo gravity high with a field of several 3D line segments. The gravity data set, given at an equidistant regular grid, was obtained from Pohánka (2001). The origin of local planar coordinates was chosen at the point of latitude 47∘57′ N, longitude 18∘00′ E, and height above sea level of 110 m, the x-coordinate being easting and y-coordinate northing. The gravity data are represented by the vertical component of the gravitation acceleration vector. Since the given area lies in lowlands and is relatively flat (variations in height are of only a few meters), no topographic correction was applied, while the data were vertically reduced to the reference altitude of 110 m using the free air gradient. The regular equally spaced grid of data was produced from the original irregularly placed data points by an inhouse interpolation method described in (ibid). The interpolation method smooths the data according to a given parameter, the so called smoothing distance. Gravity data were interpolated onto a 200 by 200meter grid. A smoothing distance of 400 m was chosen, while the mean minimum distance between the original data points was 425 m, which resulted in minimum smoothing effect. A mean value over the regular grid was subtracted from the gridded gravity. We take the residual field after subtracting the model of the regional field and after eliminating the signal of shallow sources by means of the upward and downward continuation. Three line segments are sufficient for fairly accurate approximation of the given data. The RMS of differences between the observed data and the field of the line segments is 0.57 mGal. The gravity field (without regional component) and the field of 3D line segments are presented in Fig. 3. All ends of line segments are located at depths between 6 and 9 km. It should be noted that we were seeking 21 parameters of the line segments respective to nearly 60,000 observations. Hence, this procedure is quite stable. 4. Method of local corrections for 3D data inversion In this section our approach to the 3D potential field data inversion is reviewed. We begin with the case of a restricted source body. We seek geometry of an unknown homogeneous body; the only requirement is that the sought body is assumed to be star convex relative to some point of the body. This means that there exists a point of the body such that the line segment from this point to any other point of the body is contained in the body. In this case we can introduce spherical coordinates r, θ, φ relative to the point. The body boundary can be determined by the equation r = r(θ, φ), where r(θ, φ) is a single-valued function. We introduce a Cartesian coordinate system, the origin of which coincides with the center of the spherical coordinate system. Let r = (x, y, z) be the radius-vector of an observation point. The gravitational inverse problem can be reduced to the following nonlinear equation respective to an unknown function r(θ, φ): GΔσ∫∫ K ðx; y; z; θ; φ; r ðθ; φÞÞ dθdφ ¼ U ðx; y; zÞ;

ð6Þ

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

a

475

b

20

10

15 10

5

0

km

km

5

0

−5 −10

−5

−15 −20 −20

−15

−10

−5

0

5

10

15

−10 −10

20

km 0

2

4

6

8

−5

0

5

10

km

10

12

14

16

18

20

mGal

c

d

10 mGal

20

gravity data model field

15 10

5 mGal

km

5 0 −5 3 km −10 6 km

−15

9 km -20

-15

-10

-5

0

5

10

15

−20 −20

20

−15

−10

−5

0

5

10

15

20

km

km −0.06

−0.04

−0.02

0.00

0.02

0.04

0.06

mGal Fig. 2. Approximation of synthetic data with 3D line segments. a) Synthetic gravity field of two prisms; b) and c) projections of the prisms and the obtained 3D line segments to the horizontal and vertical coordinate planes; d) residuals of approximations, RMS = 0.0188 mGal.

where G is the gravitational constant, Δσ is the density contrast, U(x, y, z) is the given field (gravitational potential, some of its derivative, or their combination). The density contrast Δσ is included in Eq. (6) as a numerical factor. It allows obtaining an equivalent set of solutions generating—with different values of the density contrast—the same gravitational field. For a fixed value of Δσ, the solution of the inverse problem is unique—for this preselected class of solutions—according to Novikov's theorem (Novikov, 1938). We propose the method of local corrections for solving nonlinear potential inverse problems similar to that given in Eq. (4). The idea is that, in computing the integrals in Eq. (6) and also in evaluating the field values on some closed surface r = ρ(θ, φ), containing the object sought within it, the very same nodes {θi, φj} are used. In each iteration an attempt is made to reduce the difference between the given and approximate field values at a fixed node solely by modifying the value of the unknown function at the same node. These considerations lead to decomposition of the inverse problem and

reduction of time expenditures to solve it, approximately by an order of magnitude. After discretization of Eq. (6) by numerical integration we obtain   cGΔσ ∑∑ K i 0 j 0 rij ¼ Ui0 j0 ; i

j

ð7Þ

where c is the weight of the cubature formula, rij = r(θi, φj), and Ui0j0 = U(θi0, φj0, ρ(θi0, φj0)), Ki0j0(rij) = K(θi0, φj0, ρ(θi0, φj0), θi, φj, rij). Our goal is to develop an iterative procedure for solving the system of nonlinear Eq. (7). Suppose that rijn are the values of the unknown function obtained at the n-th step. The idea of a local correction described above similar to the case of a contact surface, see e.g., Prutkin and Saleh (2009) leads to the fundamental equation to find the next approximation        nþ1 n n GΔσ Kij rij −Kij rij ¼ α Uij −Uij ;

ð8Þ

476

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

20

where V(r) denotes the gravitational potential of the body. For the function Ψ(r) introduced in Eq. (9) the following relation holds

10

ΨðrÞ ¼

4 8

km

12

0

4

where b = x sin θ cos φ + y sin θ sin φ + z cos θ, and a 2 = x 2 + y 2 + z 2. The integrand in Eq. (10) is algebraic relative to the unknown function and does not contain its derivatives. In the case of the core–mantle boundary (Prutkin, 2008), harmonic coefficients of the gravitational potential are given on the Earth’s surface. Based on Eq. (9), in this case we can calculate from them the harmonic coefficients of Ψ(r) and then generate its values on the Earth’s surface which we use as a closed surface containing an object with unknown geometry (the Earth's core). In the case of a local object, we benefit from the preliminary approximation with the 3D line segments. After this procedure, the parameters of segments are known. Therefore, while searching for the geometry of an object causing the same gravity effect as a given set of line segments, we can easily calculate values of Ψ(r) on any closed surface allegedly containing the unknown body. To improve the solution, one can proceed from the local corrections method to solving the integral equation for the inverse problem by Newton’s method. This method is applied rather rarely, since it ensures convergence (yet extremely rapid) only in the presence of an initial approximation close to the solution sought. In our case, such an approximation is provided by the method of local corrections. New integral equation for the magnetic inverse problem, as well as integral equations for the linearized gravitational and magnetic inverse problems can be found in Prutkin (2008). It should be noted that in the case of a contact surface we can apply the method of local corrections directly to gravity data on the physical surface. The approximation with 3D line segments is used only if we need to separate sources. The formula to find the next approximation is easier than for a restricted object. For more details we refer readers to Prutkin and Saleh (2009).

8

−10

4

−20 −20

−10

0

10

20

km −15 −11 −7

−3

1

5

9

13

17

21

8

4

25

mGal 20

10 4

km

12 12

0 20 16

8

−10

4

−20 −30

ð10Þ

12 16

−30

GΔσ 2π π r 3 ðθ; φÞsin θ dθd φ ∫0 ∫0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 r 2 ðθ; φÞ−2brðθ; φÞ þ a2

−20

0

−10

10

20

km −15 −11 −7

−3

1

5

9

13

17

21

5. Model solutions for the geological section

25

mGal Fig. 3. Approximation of the Kolarovo gravity anomaly with line segments. Top— observed data after subtracting the model regional field (a 2D harmonic function); bottom—gravitational effect of three line segments. The RMS of residuals is 0.57 mGal.

where Uijn) is the solution of the direct problem corresponding to the model rijn obtained at the n-th step. The coefficient α on the right hand side of Eq. (8) is introduced to slow down the change of the model. The integral equation of the 1st kind in Eq. (6) is ill-posed and requires regularization. Diminishing the coefficient α in Eq. (8), we can prevent highly oscillating solutions. Therefore, this factor is really similar to a regularization parameter. From Eq. (8) it is clear that time expenditures to solve the inverse problem depend dramatically on complexity of the kernel K in Eq. (6). If we use to derive Eq. (6) a representation of the field in the form of a volume integral, then the kernel function in the corresponding equation is rather complicated, and, in addition, it is a transcendental function of r (θ,φ). If we apply a representation in the form of a surface integral, then the integrand includes derivatives of the unknown function, but the differentiation of the approximately computed function r(θ,φ) is an illposed problem. We have derived new integral equations for gravitational and magnetic inverse problems which are free from both disadvantages. We introduce the following function (cf. Prutkin, 2008): 1 1 ∂V ΨðrÞ ¼ V ðrÞ− r⋅∇V ðrÞ ¼ V− r ; 2 2 ∂r

ð9Þ

After the separation of gravitational signals we transform the formal sources used as anomalous sources approximation (3D line segments) into objects of more geological meaning. For a given set of segments we obtain a restricted body or a contact surface generating the same gravity. For the Kolarovo gravity anomaly, the central line segment has substantially bigger line density than the two other. One can treat this as separation of sources and relate the field of the two lighter segments to the basement topography. We attribute the effect of the main segment to the anomalous body (a basaltic intrusion). According to results from available interpretations, we assume for the ambient medium a density of 2700 kg/m3. The density of the anomalous body equals to 3000 kg/m3. The main segment is substituted by an equivalent (relative to the external field) homogeneous body with the corresponding density contrast. The body is located below the upper boundary of the basement. Depths to the basement vary from 2 to 3 km. In Fig. 4 (top) a projection to the horizontal coordinate plane is shown, basement topography is presented by means of depth isolines. The body (with the line segment inside) is also displayed. Two other figures show projections to the vertical coordinate planes. In the upper half the solid line corresponds to given gravity, the dashed line to the model field. Depths to the center of mass, as well as to the upper and lower boundary of the body are in good agreement with the previous interpretations. By means of our methodology we can obtain alternative solutions of the geological section that are also admissible from the viewpoint of gravity data. For instance, we can attribute the whole gravity effect to a contact surface separating the underlying layer with density

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

477

are nearly flat, the observed field is dominantly caused by the shape of the source body. In this case, we obtain an intrusion with more complex geometry (see Fig. 6, top). The body is equivalent in the gravitational effect to all three line segments obtained by the line segment approximation. The last admissible solution of the geological section presented is based on the following considerations. We have applied to the given gravity the algorithm of the upward and downward continuation and observed that even if we remove sources up to the depth of 10 km, a residual signal with the amplitude of 10 mGal is still present (a half of the initial amplitude). It looks like a low-frequency isometric anomaly. We relate this component of gravity to an elevation of the dense layer below the body. If we subtract this low-frequency signal, the residual field can be approximated quite accurately (on the level of 0.75 mGal) with using only two 3D line segments. This approximation is more preferable, because two segments have lower number of parameters than three segments, so this procedure is more stable. In Fig. 6 (bottom) the last solution of the geological section is presented. It includes the contact surface (the elevation of the upper boundary of the dense layer) and the intrusion of the dense material above (it corresponds to both the 3D line segments). 6. Summary and conclusions The objective of this paper was to introduce the new methodology for 3D inversion of potential field data (with focus on gravity anomalies) in terms of signal separation both in depth and in lateral sense and the method of local corrections. To introduce the methodology we have demonstrated its capabilities on a case of the Kolarovo gravity high in the Danube Basin, southern Slovakia. The methodology belongs to the realm of direct inversion methods, making use of no forward modeling. It works with a preselected class of solutions: 20 16

18 16

20

18

10

12

km

16

14 14

20 18

0

10

12

16

8 14

18

−10

20 16 18 16

−20

18

−30 Fig. 4. Anomalous body below basement with topography. Top—projection of the body equivalent to the central 3D line segment onto the horizontal coordinate plane (the line segment is shown also), basement topography is presented by depth isolines. Middle and bottom—projections to the vertical coordinate planes. Solid line—given data, dashed line—model field.

3000 kg/m 3 from the upper layer with density 2700 kg/m 3 and estimating its depth. In this case, we do not need the approximation with 3D line segments. We invert directly the gravity data on the physical surface (without regional component). The depths of the contact surface obtained by the method of local corrections vary from 22.5 km to 7.5 km. The solution is presented in Fig. 5 in two ways: as a map of the contact surface topography (depth isolines are shown) and as a 3D surface. For visualization of 3D objects we apply Poisson Surface Reconstruction software (Kazhdan et al., 2006). The next admissible solution indicates that both the basement topography and the upper boundary of the dense layer below the body

−20

20

−10

0

10

20

km 5

7

9

11

13

15

17

19

21

23

25

km

Fig. 5. Undulations of the density contrast contact surface. Top—the obtained 3D topography of the density interface as a map of depth isolines, bottom—shaded 3D surface of the density interface.

478

I. Prutkin et al. / Journal of Applied Geophysics 75 (2011) 472–478

Fig. 6. Two additional admissible solutions. Top—intrusion with complex geometry causing the same gravity as the three line segments. Bottom—anomalous body above a dense layer with topography. The body is equivalent (by its gravitational effect) to two line segments used for the approximation.

without any linearization. This method does not make use of nonlinear minimization, which reduces the computer calculation time by an order of magnitude. 4. To find geometry of a 3D restricted object, we propose new integral equations for gravitational and magnetic inverse problems. Their integrands are algebraic relative to the unknown function and do not contain its derivatives. The use of modified potentials instead of the conventional gravitational and magnetic potentials does not introduce any additional difficulties. 5. Our methodology is successfully applied in the interpretation of the Kolarovo gravity high in Slovakia. We separated sources and inverted the gravity data by means of the herein described algorithms. This methodology provides various model solutions for the geological section: a body inside a layer with topography, an intrusion with complex geometry, an elevation of a contact surface, a body above a layer with topography. Among these admissible solutions, one can make a choice based on additional geological and geophysical constraining information. Acknowledgements

source bodies of compact star convex geometry, and 3D contact surfaces representing subsurface density interfaces. The ultimate goal was not to find the best interpretation for the Kolarovo gravity data. The goal was to demonstrate that the methodology yields a suite of admissible solutions that can be further discriminated based on additional geophysical, geologic or tectonic constraints. The quest for the best Kolarovo interpretation will be the subject of a separate paper currently under preparation. The following conclusions are drawn from our study: 1. A new algorithm is proposed to eliminate the contribution of potential field sources from the Earth’s surface to a prescribed depth h, based on applying the upward and downward continuation procedures. The solution to the 2D Dirichlet problem is used to produce a model of the regional field. Subtracting the regional field from the observed data prior to the upward continuation allows the integration of gravity data in the restricted study area. The downward continuation procedure produces the part of the field which is harmonic above the preselected depth h. The algorithm is successfully applied to separate the effects of shallow and deeper sources for the Kolarovo gravity anomaly interpretation. 2. Preliminary approximation using the 3D line segments allows separating the sources not only in depth, but also in lateral direction. Moreover, it provides reasonable mass estimates and valuable information about the shape of the source body. We divide the interpretation process into two steps: we begin with data approximation by the field of formal sources (3D line segments) with no geological meaning; in the second step we join segments and attribute their effects to a restricted body or a contact surface with the same gravity. In this latter step different solutions for the geological section can be investigated. 3. The method of local corrections is developed to recover 3D geometry of a contact surface or a restricted body. The method offers a simple and effective procedure for solving the nonlinear inverse problem

The Kolarovo gravity data interpolated to a dense regular grid at a level surface were produced by Vladimir Pohanka. Peter Vajda was supported by the Slovak Research and Development Agency under contract no. APVV-0194-10 and by Vega grant agency under project no. 2/0107/09. Miroslav Bielik was supported by Vega grant agency under project no. 1/0461/09. References Bielik, M., Fusán, O., Plančár, J., Biela, A., Túnyi, I., 1986. Some knowledge on deepseated structure of the Danube Basin (in Slovak) Geologické Práce, Správy 84, 119–134. Kazhdan, M., Bolitho, M., Hoppe, H., 2006. Poisson surface reconstruction. In: Polthier, K., Sheffer, A. (Eds.), Proceedings of the fourth Eurographics symposium on Geometry processing (SGP ’06). Eurographics Association, Aire-la-Ville, Switzerland, pp. 61–70. Kubeš, P., Bezák, V., Kucharič, L., Filo, M., Vozár, J., Konečný, V., Kohút, M., Gluch, A., 2010. Magnetic field of the Western Carpathians (Slovakia): reflections on the structure of the crust. Geologica Carpathica 61 (5), 437–447. Lavrent'ev, M.M., Romanov, V.G., Shishatskii, S.P., 1986. Ill-posed problems of mathematical physics and analysis. Translations of Math. Monographs v. 64. Amer. Math. Soc, Providence. Martyshko, P.S., Koksharov, D.E., 2005. On density determination in the layered medium by the gravity data (in Russian) Geophysical Journal 27 (4), 678–684. Martyshko, P.S., Prutkin, I.L., 2003. Technology of depth distribution of gravity field sources (in Russian) Geophysical Journal 25 (3), 159–168. Novikov, P.S., 1938. Sur le problème inverse du potentiel. Doklady Akademii Nauk SSSR 18, 165–168. Pohánka, V., 2001. Application of the harmonic inversion method to the Kolárovo gravity anomaly. Contributions to Geophysics and Geodesy 31 (4), 603–620. Prutkin, I., 2008. Gravitational and magnetic models of core–mantle boundary and their correlation. Journal of Geodynamics 45 (2–3), 146–153. Prutkin, I., Casten, U., 2009. Efficient gravity data inversion for 3D topography of a contact surface with application to the Hellenic subduction zone. Computers & Geosciences 35 (2), 225–233. Prutkin, I., Saleh, A., 2009. Gravity and magnetic data inversion for 3D topography of the Moho discontinuity in the northern Red Sea area, Egypt. Journal of Geodynamics 47 (5), 237–245. Sitárová, A., Bielik, M., Burda, M., 1984. Interpretation of the Kolárovo gravity anomaly (in Slovak) Geologické Práce, Správy 81, 171–182. Vajda, P., Bielik, M., Pohánka, V., 2002. Testing the application of the Truncation Filtering Methodology in interpreting real gravity data: the Kolárovo gravity anomaly. Contributions to Geophysics and Geodesy 32 (1), 57–66.