**New Techniques for Potential Field Data Interpretation**

Khalid S. Essa

*Cairo University/ Faculty of Science/ Geophysics Department, Giza, Egypt* 

## **1. Introduction**

20 New Achievements in Geoscience

Zouaghi, T.; Guellala, R.; Lazzez, M.; Bédir, M.; Ben Youssef, M.; Inoubli, M.H. & Zargouni,

Zouari, H. ; Turki, M. M. & Delteil J. (1990). Nouvelles données sur l'évolution tectonique de la chaîne de Gafsa. *Bulletin de la Société Géologique de France 8*, pp. 621-628. Zouari, H. (1995). Evolution géodynamique de l'Atlas centro-méridional de la Tunisie.

Zouhri, L.; Gorini, C.; Lamouroux, C.; Vachard, D.& Dakki M. (2003). Interprétation

Zouhri, L.; Gorini, C.; Mania, J; Deffontaines, B. & Zerouali, A. (2004). Spatial distribution of

in the Rharb basin, Morocco: *Hydrological Science Journal*, 49, 387–398.

Austria.

*Doct. es-Sciences, Univ. Tunis II,* 278 p.

sismique réflexion. *C.R. Acad. Sci* . pp 319-326.

F. (2011). *The Chotts Fold Belt of Southern Tunisia, North African Margin: Structural Pattern, Evolution, and Regional Geodynamic Implications. New Frontiers in Tectonic Research - At the Midst of Plate Convergence. ISBN: 978-953-307-594-5. Intech* Vienna,

Stratigraphie, analyses géométrique, cinématique et tectono-sédimentaire. *Thèse* 

hydrogéologique de l'aquifère des bassins Sud-rifains (Maroc): apport de la

resistivity in the hydrogeological systems, and identification of the catchment area

Many of the geological structures in mineral and petroleum exploration can be classified into four categories: spheres, cylinders, dikes and geological contacts. These four simple geometric forms are convenient approximations to common geological structures often encountered in the interpretation of magnetic and gravity data. Two quantitative categories are usually adapted for the interpretation of magnetic and gravity anomalies in order to determine basically the depth, the shape of the buried structure and the other model parameters. First, two parameters are considered the most important problem in exploration geophysics. The first category includes 2D and 3D continuous modelling and inversion methods (Mohan et al., 1982; Nabighian, 1984; Chai & Hinze, 1988; Martin-Atienza & Garcia-Abdeslem 1999; Zhang et al., 2001; Keating & Pilkington, 2004; Numes et al., 2008; Martins et al., 2010; Silva et al., 2010). The solutions of these methods require magnetic susceptibility and/or remanent magnetization for magnetic interpretation and density information for gravity interpretation a as part of the input, along with the same depth information obtained from geological and/or other geophysical data. Thus, the resulting model can vary widely depending on these factors because the magnetic and gravity inverse problem are ill-posed and are, therefore unstable and nonunique (Zhdanov, 2002; Tarantola, 2005). Standard techniques for obtaining a stable solution of an ill-posed inverse problem include those based on regularization methods (Tikhonov & Arsenin, 1977). The second category is fixed simple geometry methods, in which the spheres, cylinders, dikes, and geological contacts models estimate the depth and the shape of the buried structures from residuals and/or observed data. The models may be shifted from geological reality, but they are usually sufficient to determine whether the form and the magnitude of both calculated magnetic and gravity effects are close enough to those observed to make the geological interpretation reasonable. Several methods have been developed for the second category to interpret magnetic and gravity data using a fixed simple geometry (Nettleton, 1976; Stanley, 1977; Atchuta Rao & Ram Babu, 1980; Prakasa Rao & Subrahmanyan, 1988; Abdelrahman & El-Araby, 1993; Zhang et al., 2000; Salem et al., 2004; Abdelrahman & Essa, 2005; Essa, 2007; Essa, 2011). The drawbacks of these methods are that they are highly subjective, need a priori information about the shape (shape factor) of the anomalous body and use few characteristic points when inverting for the remaining parameters. New algorithms for magnetic and gravity inversion have been developed that successively determines the depth

New Techniques for Potential Field Data Interpretation 23

Table 1. Definition of a, b, c, m, n, p, r, and q values were shown in equation (1). F.H.D. and

Substituting equation (4) into equation (3), we obtain the following nonlinear equation in z

*NT z az bx ax T N N z x T z az bN Tx z*

(0) ( ) ( )( ) (0) ( ) ( ,) , ( ) *q r r r q q r i i i*

*i*

*i i*

() ( ) ( ,) ,

*aN x z* 

2 2 2 2 2 2 2 2 2 2 2 2

2

*z Lx Tx z* (6)

(5)

S.H.D. are the first and the second horizontal derivatives of the magnetic anomaly,

*N*

*i*

*i q*

The unknown depth z in equation (5) can be obtained by minimizing

where L (xi) denotes the observed magnetic anomaly at xi.

Fig. 1. Various simple geometrical structures diagram.

respectively.

(z) and the other associated model parameters of the buried structure. The inverse problem of the depth estimation from the observed magnetic and gravity data has been transformed into a nonlinear equation of the form f(z) = 0. This equation is then solved for depth by minimizing an objective functional in the least-squares sense. Using the estimated depth, the other model parameters are computed from the measured magnetic and gravity data, respectively. The procedures are applied to synthetic data with and without noise for magnetic and gravity data. These procedures are also tested to complicated regional and interference from neighbouring magnetic rocks. Finally, these techniques are also successfully applied to real data sets for mineral and petroleum exploration, and it is found that the estimated depths and the associated model parameters are in good agreement with the actual values.

#### **2. Formulation of the problems**

#### **2.1 A least-squares minimization approach to depth determination from residual magnetic anomaly**

Following Gay (1963), Prakasa Rao et al. (1986), and Prakasa Rao & Subrahmanyan (1988), the magnetic anomaly expression produced by most geologic structures with center located at xi= 0 can be represented by the following equation:

$$\mathbf{T(x\_i, x, \theta) = K} \frac{\left(\mathbf{a}\mathbf{z}^{\mathbf{2r}} + \mathbf{b}\mathbf{x\_i^{\mathbf{2}}}\right) \left(\sin\theta\right)^{\mathbf{m}} \left(\cos\theta\right)^{\mathbf{n}} + \mathbf{c}\mathbf{x\_i}\mathbf{z}^{\mathbf{P}} \left(\sin\theta\right)^{\mathbf{n}} \left(\cos\theta\right)^{\mathbf{m}}}{\left(\mathbf{x\_i^{2}} + \mathbf{z}^{2}\right)^{\mathbf{q}}}, \mathbf{i = 1, 2, 3, \dots \text{L}}, \mathbf{i = 1, 2, 3, \dots \text{L}} \qquad \text{(1)}$$

The geometries are shown in Figure 1. In equation (1), z is the depth, K is the amplitude coefficient (effective magnetization intensity), θ is the index parameter (effective magnetization inclination), xi is the horizontal coordinate position, and q is the shape factor. Values for a, b, c, m, n, r, p, and q, are given in Table 1.

At the origin (xi =0), equation (1) gives the following relationship:

$$K = \frac{T(0)z^{2q-2r}}{a(\sin \theta)^m (\cos \theta)^n},\tag{2}$$

where T(0) is the anomaly value at the origin (Fig. 2).

Using equation (2), equation (1) can be rewritten as

$$T(X\_i, Z, \theta) = \frac{T(0)z^{2q-2r}}{a} \left[ \frac{(az^{2r} + bx\_i^2) + cx\_i z^p (\tan \theta)^{n-m}}{(x\_i^2 + z^2)^q} \right],\tag{3}$$

For all shapes (function of q), equation (3) gives the following relationship at xi=N

$$(\tan \theta)^{n-m} = \frac{aT(N)(N^2 + z^2)^q - T(0)z^{2q-2r}(az^{2r} + bN^2)}{cNz^pT(0)z^{2q-2r}},\tag{4}$$

where T (N) is the anomaly value at xi = N (Fig. 2).

(z) and the other associated model parameters of the buried structure. The inverse problem of the depth estimation from the observed magnetic and gravity data has been transformed into a nonlinear equation of the form f(z) = 0. This equation is then solved for depth by minimizing an objective functional in the least-squares sense. Using the estimated depth, the other model parameters are computed from the measured magnetic and gravity data, respectively. The procedures are applied to synthetic data with and without noise for magnetic and gravity data. These procedures are also tested to complicated regional and interference from neighbouring magnetic rocks. Finally, these techniques are also successfully applied to real data sets for mineral and petroleum exploration, and it is found that the estimated depths and the associated model parameters are in good agreement with

**2.1 A least-squares minimization approach to depth determination from residual** 

Following Gay (1963), Prakasa Rao et al. (1986), and Prakasa Rao & Subrahmanyan (1988), the magnetic anomaly expression produced by most geologic structures with center located

 

The geometries are shown in Figure 1. In equation (1), z is the depth, K is the amplitude coefficient (effective magnetization intensity), θ is the index parameter (effective magnetization inclination), xi is the horizontal coordinate position, and q is the shape factor.

> 2 2 (0) , (sin ) (cos ) *q r m n*

 

2 2

*q r r p n m i i*

*i*

2 2 2 2 2 2 2 2

(4)

*p qr aT N N z T z az bN cNz T z*

(2)

2 2 2 2

*i q*

For all shapes (function of q), equation (3) gives the following relationship at xi=N

*T z az bx cx z TX Z*

(0) ( ) (tan ) ( ,,) , ( )

( )( ) (0) ( ) (tan ) , (0) *q qr <sup>r</sup> n m*

*a x z*

*T z <sup>K</sup> a*

(1)

(3)

 2r 2 mn nm <sup>p</sup> az bx sin<sup>θ</sup> cos<sup>θ</sup> cx z sin<sup>θ</sup> cos<sup>θ</sup> i i T(x ,z ,<sup>θ</sup> ) K ,i 1,2,3,...L <sup>i</sup> <sup>q</sup> <sup>2</sup> <sup>2</sup> x z <sup>i</sup>

the actual values.

**magnetic anomaly** 

**2. Formulation of the problems** 

at xi= 0 can be represented by the following equation:

Values for a, b, c, m, n, r, p, and q, are given in Table 1.

where T(0) is the anomaly value at the origin (Fig. 2). Using equation (2), equation (1) can be rewritten as

where T (N) is the anomaly value at xi = N (Fig. 2).

At the origin (xi =0), equation (1) gives the following relationship:

Fig. 1. Various simple geometrical structures diagram.


Table 1. Definition of a, b, c, m, n, p, r, and q values were shown in equation (1). F.H.D. and S.H.D. are the first and the second horizontal derivatives of the magnetic anomaly, respectively.

Substituting equation (4) into equation (3), we obtain the following nonlinear equation in z

$$T(\mathbf{x}\_i, \mathbf{z}) = \frac{NT(0)z^{2q-2r}(az^{2r} + b\mathbf{x}\_i^2) + az\_i T(\mathbf{N})(\mathbf{N}^2 + z^2)^q - \mathbf{x}\_i T(0)z^{2q-2r}(az^{2r} + b\mathbf{N}^2)}{aN(\mathbf{x}\_i^2 + z^2)^q},\tag{5}$$

The unknown depth z in equation (5) can be obtained by minimizing

$$\boldsymbol{\Psi}(\mathbf{z}) = \sum\_{i}^{N} \left[ L(\mathbf{x}\_{i}) - T(\mathbf{x}\_{i}, \mathbf{z}) \right]^{2},\tag{6}$$

where L (xi) denotes the observed magnetic anomaly at xi.

New Techniques for Potential Field Data Interpretation 25

should be restricted to as few parameters as is consistent with obtaining useful results. This is why we propose a solution for only one unknown, z. Once z is known, the effective magnetization inclination, θ, can be determined from equation (4). Finally, knowing θ, the effective magnetization intensity, K, can be determined from equation (2). Then, we measure the goodness of fit between the observed and the computed magnetic data for each N value. The simplest way to compare two magnetic profiles is to compute the root-mean-square (RMS) of the differences between the observed and the fitted anomalies. The model parameters which give the least root-mean square error are the best. To this point we have assumed knowledge of the axes of the magnetic profile so that T(0) can be found. T(0) is determined using methods described by Stanley (1977). As illustrated in Figure 2, the line M-m intersects the anomaly profile at xi = 0. The base line of the anomaly profile lies a distance M-T(0) above the minimum. A semi-automated interpretation scheme based on the

above equations for analyzing field data is illustrated in Figure 3.

Fig. 3. Generalized scheme for semi-automatic depth, index parameter and amplitude

Synthetic examples of spheres and horizontal cylinders buried at different depths (profile length = 40 units, K = 100 units, θ = 60o, sampling interval = 1unit) were interpreted using the present method [equations (7), (4), and (2)] to determine depth, index parameter, and amplitude coefficient, respectively. In each case, the starting depth was 5 units. In all cases

coefficient estimations.

**2.1.1 Synthetic example** 

Fig. 2. A typical magnetic anomaly profile over a thin dike. The anomaly value at the origin T(0) and the anomaly value T(N) where N is taken to be 1arbitrary unit in this case, the position of the maximum value (M), and the minimum value (m) are illustrated.

Minimization of ψ(z) in the least - squares sense, i.e., (d/dz) ψ(z)=0 leads to the following equation:

$$f(\mathbf{z}) = \sum\_{i}^{N} \left[ L(\mathbf{x}\_{i}) - T(\mathbf{x}\_{i}, \mathbf{z}) \right] T^\*(\mathbf{x}\_{i}, \mathbf{z}) = \mathbf{0},\tag{7}$$

where

$$T^\*(\mathbf{x}\_{i\prime}z) = (\bigvee\_{d\prime z} \mathbf{\color{red}{T}(x\_{i\prime}z)} $$

Equation (7) can be solved for z using the standard methods for solving nonlinear equations. Here, it is solved by an iteration method (Press et al., 1986). The iteration form of equation (7) is given as

$$z\_f = f(z\_j) \quad , \tag{8}$$

where zj is the initial depth and zf is the revised depth. zf will be used as zj for next iteration. The iteration stops when **|zf - zj |< e**, where e is a small predetermined real number close to zero. The source body depth is determined by solving one nonlinear equation in z. Any initial guess for z works well because there is only one minimum. The experience with the minimization technique for two or more unknowns is that it always produced good results for synthetic data with or without random noise. In the case of field data, good results may only be obtained when using very good initial guesses on the model parameters. The optimization problem for the depth parameter is highly nonlinear, increasing the number of parameters to be solved simultaneously also increases the dimensionality of the energy surface, thereby greatly increasing the probability of the optimization stalling in a local minimum on that surface. Thus common sense dictates that the nonlinear optimization

Fig. 2. A typical magnetic anomaly profile over a thin dike. The anomaly value at the origin T(0) and the anomaly value T(N) where N is taken to be 1arbitrary unit in this case, the position of the maximum value (M), and the minimum value (m) are illustrated.

Minimization of ψ(z) in the least - squares sense, i.e., (d/dz) ψ(z)=0 leads to the following

*N*

*i*

\*

) ( ) ( , ) ( , ) 0,

( ,) ( )( ,) *i i T x z Tx z <sup>d</sup>*

Equation (7) can be solved for z using the standard methods for solving nonlinear equations. Here, it is solved by an iteration method (Press et al., 1986). The iteration form of equation

where zj is the initial depth and zf is the revised depth. zf will be used as zj for next iteration. The iteration stops when **|zf - zj |< e**, where e is a small predetermined real number close to zero. The source body depth is determined by solving one nonlinear equation in z. Any initial guess for z works well because there is only one minimum. The experience with the minimization technique for two or more unknowns is that it always produced good results for synthetic data with or without random noise. In the case of field data, good results may only be obtained when using very good initial guesses on the model parameters. The optimization problem for the depth parameter is highly nonlinear, increasing the number of parameters to be solved simultaneously also increases the dimensionality of the energy surface, thereby greatly increasing the probability of the optimization stalling in a local minimum on that surface. Thus common sense dictates that the nonlinear optimization

*ii i*

*dz*

*f(z L x T x z T x z* (7)

( ), *f j z fz* (8)

equation:

where

(7) is given as

should be restricted to as few parameters as is consistent with obtaining useful results. This is why we propose a solution for only one unknown, z. Once z is known, the effective magnetization inclination, θ, can be determined from equation (4). Finally, knowing θ, the effective magnetization intensity, K, can be determined from equation (2). Then, we measure the goodness of fit between the observed and the computed magnetic data for each N value. The simplest way to compare two magnetic profiles is to compute the root-mean-square (RMS) of the differences between the observed and the fitted anomalies. The model parameters which give the least root-mean square error are the best. To this point we have assumed knowledge of the axes of the magnetic profile so that T(0) can be found. T(0) is determined using methods described by Stanley (1977). As illustrated in Figure 2, the line M-m intersects the anomaly profile at xi = 0. The base line of the anomaly profile lies a distance M-T(0) above the minimum. A semi-automated interpretation scheme based on the above equations for analyzing field data is illustrated in Figure 3.

Fig. 3. Generalized scheme for semi-automatic depth, index parameter and amplitude coefficient estimations.

#### **2.1.1 Synthetic example**

Synthetic examples of spheres and horizontal cylinders buried at different depths (profile length = 40 units, K = 100 units, θ = 60o, sampling interval = 1unit) were interpreted using the present method [equations (7), (4), and (2)] to determine depth, index parameter, and amplitude coefficient, respectively. In each case, the starting depth was 5 units. In all cases

New Techniques for Potential Field Data Interpretation 27

Fig. 5. Total magnetic anomaly (above) over an outcropping diabase dike (below), Pishabo Lake, Ontario, Canada (McGrath & Hood, 1970). The base line and zero crossing shown are

**2.2 Three least-squares minimization approach to model parameters estimation from** 

From equation (1), the magnetic anomaly expression of a thin dike which extends to infinity

2 2 sin cos ( ,,) , 1,2,3,....., *<sup>i</sup>*

where z is the depth to the top of the body, xi is a distance point along x-axis where the observed anomaly is located, K is the amplitude coefficient, and θ is the index parameter.

(9)

*i x z T x z zK i N x z* 

determined using Stanley's method (1977).

Table 2. Numerical results of the Pishabo field example.

in both strike direction and down dip (2-D) is given by:

*i*

**magnetic anomaly due to thin dikes** 

examined, the exact values of z, θ, and K were obtained. However, in studying the error response of the least-squares method, synthetic examples contaminated with 5% random errors were considered. Following the interpretation scheme, values of the most appropriate model parameters (z, θ, K) were computed and percentage of error in model parameters was plotted against the model depth for comparison (Fig. 4). We verified numerically that the depth obtained is within 3% for horizontal cylinders and 4% for spheres. The index parameter obtained is within 5.5% whereas the amplitude coefficient is within 8% (Fig. 4). Good results are obtained by using the present algorithm, particularly for depth estimation, which is a primary concern in magnetic prospecting and other geophysical work.

Fig. 4. Model parameters error response estimates. Abscissa: model depth. Ordinate: percent error in model parameters.

#### **2.1.2 Field example**

Figure 5 shows a total magnetic anomaly above an olivine diabase dike, Pishabo Lake, Ontario (McGrath & Hood, 1970). The geological cross section of the dike is shown beneath the anomaly. The depth to the outcropping dike (sensor height) is 304 m (Fig. 5). The anomaly profile was digitized at an interval of 100 m. Equations (7), (4), and (2) were used to determine depth, index parameter and amplitude coefficient, respectively, using all possible cases on N values. The starting depth used in this field example was 50 m. Then, we computed the rootmean-square of the differences between the observed and the fitted anomalies. The best fit model parameters are: z = 306 m, θ = 38o and K = 1422 nT\*Z (Table 2). McGrath & Hood (1970) applied a computer curve-matching method to the same magnetic data employing a leastsquares method and obtained a depth of 301 m. Moreover, Figure 5 shows that there is a lack of ''thinness'' because the thickness of the dike is only slightly less than the depth. However, because many of the characteristics of the thick dike anomalies are similar to those for thin dike (Gay, 1963), our method can be applied not only to magnetic anomalies due to thin dikes, but also to those due to thick dikes to obtain reliable depth estimates.

examined, the exact values of z, θ, and K were obtained. However, in studying the error response of the least-squares method, synthetic examples contaminated with 5% random errors were considered. Following the interpretation scheme, values of the most appropriate model parameters (z, θ, K) were computed and percentage of error in model parameters was plotted against the model depth for comparison (Fig. 4). We verified numerically that the depth obtained is within 3% for horizontal cylinders and 4% for spheres. The index parameter obtained is within 5.5% whereas the amplitude coefficient is within 8% (Fig. 4). Good results are obtained by using the present algorithm, particularly for depth estimation,

Fig. 4. Model parameters error response estimates. Abscissa: model depth. Ordinate: percent

Figure 5 shows a total magnetic anomaly above an olivine diabase dike, Pishabo Lake, Ontario (McGrath & Hood, 1970). The geological cross section of the dike is shown beneath the anomaly. The depth to the outcropping dike (sensor height) is 304 m (Fig. 5). The anomaly profile was digitized at an interval of 100 m. Equations (7), (4), and (2) were used to determine depth, index parameter and amplitude coefficient, respectively, using all possible cases on N values. The starting depth used in this field example was 50 m. Then, we computed the rootmean-square of the differences between the observed and the fitted anomalies. The best fit model parameters are: z = 306 m, θ = 38o and K = 1422 nT\*Z (Table 2). McGrath & Hood (1970) applied a computer curve-matching method to the same magnetic data employing a leastsquares method and obtained a depth of 301 m. Moreover, Figure 5 shows that there is a lack of ''thinness'' because the thickness of the dike is only slightly less than the depth. However, because many of the characteristics of the thick dike anomalies are similar to those for thin dike (Gay, 1963), our method can be applied not only to magnetic anomalies due to thin dikes,

but also to those due to thick dikes to obtain reliable depth estimates.

error in model parameters.

**2.1.2 Field example** 

which is a primary concern in magnetic prospecting and other geophysical work.

Fig. 5. Total magnetic anomaly (above) over an outcropping diabase dike (below), Pishabo Lake, Ontario, Canada (McGrath & Hood, 1970). The base line and zero crossing shown are determined using Stanley's method (1977).


Table 2. Numerical results of the Pishabo field example.

#### **2.2 Three least-squares minimization approach to model parameters estimation from magnetic anomaly due to thin dikes**

From equation (1), the magnetic anomaly expression of a thin dike which extends to infinity in both strike direction and down dip (2-D) is given by:

$$T(\mathbf{x}\_i, z, \theta) = zK \frac{\mathbf{x}\_i \sin \theta + z \cos \theta}{\mathbf{x}\_i^2 + z^2}, \text{ i = 1, 2, 3, ..., N} \tag{9}$$

where z is the depth to the top of the body, xi is a distance point along x-axis where the observed anomaly is located, K is the amplitude coefficient, and θ is the index parameter.

New Techniques for Potential Field Data Interpretation 29

2 2 2 cos (0) . *Ks <sup>R</sup> s z*

( ) (0) 2 tan 2 ( )tan ( )tan ( ,, ,) . <sup>2</sup> () ()

2 2 ( ) 2 3 tan , (0) <sup>4</sup> *<sup>i</sup>*

2 2 ( ) 2 3 tan , . (0) <sup>4</sup> *<sup>i</sup>*

> 2 2 (4 ) tan , <sup>6</sup> *Fs z zs*

() ( ) , (0) *Rs R s <sup>F</sup> R* 

2 2 22 2 22 2 3 2 2 2 2

*s xz xs z*

*i i*

( ) ( ) (0) ( , , ) ,

*z Lx R W x zs*

*i i*

*i*

which is computed directly from observed anomaly. In all cases, the maximum value of s is less than or equal (N-1)/4, where N is the number of the data points along the anomaly profile whose centre is located at xi = 0. Using equation (19), equation (16) can be written as

( ) 2 (4 ) 12 ( ) (4 ) 6 ( ,,) <sup>12</sup> ( )

The unknown depth (z) in equation (20) can be obtained by minimizing

1

*i*

*N*

*s z xF s z z s x sF s z z s Wx zs*

*i i*

*zs z R x z xs z xs z Rx z s*

2 2

2 2

*R s z s zs*

*R s z* 

*R s z s zs*

*R s z*

2 22 2 2 2 2

*s x z xs z xs z*

(16)

*ii i*

*x s*

*x s*

(19)

( , , ) (0) ( , , ), *Rx zs R W x zs i i* (20)

22 2 2 2

( ) (4 ) 6 . ( )

*x sF s z zs xs z*

*i*

2

(21)

*ii i*

(15)

(17)

(18)

The specific value R(xi, z, θ, s) at xi = 0 is given by

Using equation (15), equation (14) can now be written as

Using equation (16), we can obtain the following equations at xi = s

2 2

*i*

And

where

where

*i*

From equations (17) and (18), we obtain

The values of K and θ for the anomalies in total, vertical, and horizontal fields are given in Table 3. The index parameter θ is related to the effective direction of the resultant magnetization, the dip and strike of the dike.


Table 3. Characteristic index parameter θ and amplitude coefficient K in vertical (ΔV), horizontal (ΔH), and total (ΔT) magnetic anomalies due to thin dikes where t is the thickness, k is the magnetic susceptibility contrast, d is the dip angle, To and Ío are the values of effective total intensity and effective inclination of magnetic polarization, and α is the strike (after Gay, 1963).

Following Griffin (1949), the moving average residual magnetic anomaly (R) is defined as

$$R = T(0) - \overline{T}(\text{s}) \tag{10}$$

where T(0) is the magnetic anomaly value at a given point on the magnetic anomaly map, and

$$\overline{T}(\mathbf{s}) = \frac{1}{2\pi} \int\_0^{2\pi} T(\mathbf{s}, \boldsymbol{\beta}) \, d\boldsymbol{\beta} \tag{11}$$

is the average magnetic value at the radial distance s from the point where the magnetic value is T(0). The geometrical interpretation of equation (11) is that the average value, *T s*( ) , is obtained by going around a circle of radius s, having T(0) as its centre, and forming the sum of the products of T(s, β) and dβ for an infinite number of infinitesimally small dβ's. The result of the above is then divided by 2π. Since an integral form of *T s*( ) is not know, Griffin (1949) has adopted a simple numerical method to obtain the average magnetic anomaly value. The method is to form the arithmetical average of a finite number of points about the circumference of a circle of radius s, or

$$\overline{T}(\mathbf{s}) = \frac{\left[T\_1(\mathbf{s}) + T\_2(\mathbf{s}) + T\_3(\mathbf{s}) + \dots \dots + T\_n(\mathbf{s})\right]}{n}.\tag{12}$$

for n points.

Along a profile, the moving average residual magnetic anomaly is then defined as

$$R(\mathbf{x}\_{i'}, \mathbf{s}) = \frac{2T(\mathbf{x}\_i) - T(\mathbf{x}\_i + \mathbf{s}) - T(\mathbf{x}\_i - \mathbf{s})}{2},\tag{13}$$

where s = 1, 2, 3,…., M sample spacing units and is defined here as the window length or graticule spacing.

Applying equation (13) to equation (9), the moving average residual magnetic anomaly due to a thin dike is obtained as

$$R(\mathbf{x}\_i, z, \theta, \mathbf{s}) = \frac{zK}{2} \left\{ \frac{2\mathbf{x}\_i \sin \theta + 2z \cos \theta}{\mathbf{x}\_i^2 + z^2} - \frac{(\mathbf{x}\_i + \mathbf{s}) \sin \theta + z \cos \theta}{(\mathbf{x}\_i + \mathbf{s})^2 + z^2} - \frac{(\mathbf{x}\_i - \mathbf{s}) \sin \theta + z \cos \theta}{(\mathbf{x}\_i - \mathbf{s})^2 + z^2} \right\}. \tag{14}$$

The specific value R(xi, z, θ, s) at xi = 0 is given by

$$R(0) = \frac{Ks^2 \cos \theta}{s^2 + z^2}. \tag{15}$$

Using equation (15), equation (14) can now be written as

$$R(\mathbf{x}\_i, z\_i, \theta, s) = \frac{z(s^2 + z^2)R(0)}{2s^2} \left\{ \frac{2\mathbf{x}\_i \tan \theta + 2z}{\mathbf{x}\_i^2 + z^2} - \frac{(\mathbf{x}\_i + \mathbf{s})\tan \theta + z}{\left(\mathbf{x}\_i + \mathbf{s}\right)^2 + z^2} - \frac{(\mathbf{x}\_i - \mathbf{s})\tan \theta + z}{\left(\mathbf{x}\_i - \mathbf{s}\right)^2 + z^2} \right\}.\tag{16}$$

Using equation (16), we can obtain the following equations at xi = s

$$\frac{R(s)}{R(0)} = \frac{z^2 - 2s^2 + 3zs\tan\theta}{4s^2 + z^2}, x\_i = s \tag{17}$$

And

28 New Achievements in Geoscience

The values of K and θ for the anomalies in total, vertical, and horizontal fields are given in Table 3. The index parameter θ is related to the effective direction of the resultant

Table 3. Characteristic index parameter θ and amplitude coefficient K in vertical (ΔV), horizontal (ΔH), and total (ΔT) magnetic anomalies due to thin dikes where t is the thickness, k is the magnetic susceptibility contrast, d is the dip angle, To and Ío are the values of effective total intensity and effective inclination of magnetic polarization, and α is

Following Griffin (1949), the moving average residual magnetic anomaly (R) is defined as

where T(0) is the magnetic anomaly value at a given point on the magnetic anomaly map, and 2

0 <sup>1</sup> () (, ) <sup>2</sup> *Ts Ts d* 

is the average magnetic value at the radial distance s from the point where the magnetic value is T(0). The geometrical interpretation of equation (11) is that the average value, *T s*( ) , is obtained by going around a circle of radius s, having T(0) as its centre, and forming the sum of the products of T(s, β) and dβ for an infinite number of infinitesimally small dβ's. The result of the above is then divided by 2π. Since an integral form of *T s*( ) is not know, Griffin (1949) has adopted a simple numerical method to obtain the average magnetic anomaly value. The method is to form the arithmetical average of a finite number of points

> <sup>123</sup> ( ) ( ) ( ) ....... ( ) ( ) . *Ts Ts Ts Ts <sup>n</sup> T s n*

2( ) ( ) ( ) ( ,) , <sup>2</sup> *ii i*

where s = 1, 2, 3,…., M sample spacing units and is defined here as the window length or

Applying equation (13) to equation (9), the moving average residual magnetic anomaly due

2 sin 2 cos ( )sin cos ( )sin cos ( ,, ,) . <sup>2</sup> () () *ii i*

*zK x z xs z xs z Rx z s*

*ii i*

Along a profile, the moving average residual magnetic anomaly is then defined as

*i*

  

*R T Ts* (0) ( ), (10)

(11)

(12)

*Tx Tx s Tx s Rx s* (13)

2 2 2 2 2 2

(14)

*x z xs z xs z*

magnetization, the dip and strike of the dike.

about the circumference of a circle of radius s, or

the strike (after Gay, 1963).

for n points.

graticule spacing.

*i*

to a thin dike is obtained as

$$\frac{R(-s)}{R(0)} = \frac{z^2 - 2s^2 - 3zs\tan\theta}{4s^2 + z^2}, \text{x}\_i = -s. \tag{18}$$

From equations (17) and (18), we obtain

$$\tan \theta = \frac{F(4s^2 + z^2)}{6zs},\tag{19}$$

where

$$F = \frac{R(s) - R(-s)}{R(0)}\,\prime$$

which is computed directly from observed anomaly. In all cases, the maximum value of s is less than or equal (N-1)/4, where N is the number of the data points along the anomaly profile whose centre is located at xi = 0. Using equation (19), equation (16) can be written as

$$R(\mathbf{x}\_{i\prime}, \mathbf{z}\_{\prime}\mathbf{s}) = R(\mathbf{0})W(\mathbf{x}\_{i\prime}, \mathbf{z}\_{\prime}\mathbf{s})\_{\prime} \tag{20}$$

where

$$\begin{split} W(\mathbf{x}\_{i},\mathbf{z},\mathbf{s}) = \frac{\left(\mathbf{s}^{2} + \mathbf{z}^{2}\right)}{12\mathbf{s}^{3}} \Bigg\{ \frac{2\mathbf{x}\_{i}F(4\mathbf{s}^{2} + \mathbf{z}^{2}) + 12\mathbf{z}^{2}\mathbf{s}}{\mathbf{x}\_{i}^{2} + \mathbf{z}^{2}} - \frac{\left(\mathbf{x}\_{i} + \mathbf{s}\right)F(4\mathbf{s}^{2} + \mathbf{z}^{2}) + 6\mathbf{z}^{2}\mathbf{s}}{\left(\mathbf{x}\_{i} + \mathbf{s}\right)^{2} + \mathbf{z}^{2}} \\ & - \frac{\left(\mathbf{x}\_{i} - \mathbf{s}\right)F(4\mathbf{s}^{2} + \mathbf{z}^{2}) + 6\mathbf{z}^{2}\mathbf{s}}{\left(\mathbf{x}\_{i} - \mathbf{s}\right)^{2} + \mathbf{z}^{2}} \Bigg\}. \end{split}$$

The unknown depth (z) in equation (20) can be obtained by minimizing

$$\varphi(\mathbf{z}) = \sum\_{i=1}^{N} \left[ L(\mathbf{x}\_i) - R(\mathbf{0}) W(\mathbf{x}\_i, \mathbf{z}, \mathbf{s}) \right]^2,\tag{21}$$

New Techniques for Potential Field Data Interpretation 31

Equation (25) can be solved using a simple iterative method (Demodivch & Maron, 1973). Substituting the computed depth (zc) and index parameter (θc) in equation (14) as fixed

2 sin 2 cos ( )sin cos ( )sin cos ( ,) . <sup>2</sup> ( ) ( ) *c i c c c i cc c i cc c*

Finally, applying the least-squares method to equation (26), the unknown amplitude

*z x z xs z xs z Dx s*

1

*<sup>i</sup> <sup>c</sup> <sup>N</sup>*

 

*K*

*N*

1

Theoretically, one value of s is enough to determine the model parameters (z, θ, K) from equations (22), (25), and (27), respectively. In practice, more than one value of s is desirable because of the presence of noise, interference from neighbouring magnetic rocks, and complicated regional in data. Moreover, the advantage of this method over continuous modelling methods is that it requires neither magnetic susceptibility contrast nor depth information, and it can be applied if little or no factual information other than the magnetic data is available. The method may be automated to solve for a number of anomalies along a magnetic profile when the origin of each dike is known from other geophysical and/or

The composite vertical magnetic anomaly in nanoteslas of Figure 6 consisting of the combined effect of an intermediate structure of interest (thin dike with K = 200 nT, z = 5 km, and θ = 40o), an interference neighbouring magnetic rocks (horizontal cylinder with K = 40 nT, z = 2 km, and θ = 30o), and a deep-seated structure (sphere with K = 80000 nT, z = 15 km,

2 22

2

*i*

*x*

(450 ( 60) )sin 10 45( 60)cos10 <sup>80000</sup> (( 60) 15

*x x*

*i i*

( )

*Horizontalcylinder verticalcomponent*

2

.

(28)

<sup>5</sup> <sup>2</sup> <sup>2</sup>

)

( )

2

(4 ( 10) )cos30 4( 10)sin 30 ( ) 40 (( 10) 2 )

*i*

*x*

*i i*

2 2

5

*Thindike verticalcomponent*

*Sphere verticalcomponent* ( )

sin 40 5cos 40 <sup>200</sup>

*x x T x*

*i*

*x*

*i*

*x*

*i*

 

2 2 2 2 2 2

*i c i c i c*

*x z xs z xs z* 

*Dx s*

( ,)

*i*

( ) ( ,)

*i i*

*Lx Dx s*

2

.

( , ) ( , ), *R x s KD x s i i* (26)

(27)

parameters, we obtain

coefficient (Kc) can be determined from

*i*

geological data.

**2.2.1 Synthetic example** 

and θ = 10o) was computed by the following expression:

*i*

where

where L(xi) denotes the moving average residual magnetic anomalies at xi. Setting the derivative of φ(z) to zero with respect to z leads to

$$f(\mathbf{z}) = \sum\_{i=1}^{N} \left[ L(\mathbf{x}\_i) - R(0) \mathcal{W}(\mathbf{x}\_i, \mathbf{z}, \mathbf{s}) \right] \mathbb{W}^\*(\mathbf{x}\_i \mathbf{z}, \mathbf{s}) = \mathbf{0}. \tag{22}$$

where W\*(xi, z, s) is evaluated by analytically differentiating W(xi, z, s) and is given by the following equation

$$\mathcal{W}^{\*}\left(\mathbf{x}\_{i},z,s\right) = \left(\mathbf{s}^{2} + z^{2}\right) \frac{\frac{\left(\mathbf{x}\_{i}^{2} + z^{2}\right)\left(4\mathbf{x}\_{i}Fz + 24zs\right) - 2z\left(2\mathbf{x}\_{i}F\left(4\mathbf{s}^{2} + z^{2}\right) + 12z^{2}s\right)}{\left(\mathbf{x}\_{i}^{2} + z^{2}\right)^{2}} - $$

$$\frac{\left(\left(\mathbf{x}\_{i} + s\right)^{2} + z^{2}\right)\left(2\left(\mathbf{x}\_{i} + s\right)Fz + 12zs\right) - 2z\left(\left(\mathbf{x}\_{i} + s\right)F\left(4s^{2} + z^{2}\right) + 6z^{2}s\right)}{\left(\left(\mathbf{x}\_{i} + s\right)^{2} + z^{2}\right)^{2}} - 1}{\left(\left(\mathbf{x}\_{i} + s\right)^{2} + z^{2}\right)^{2}} - \frac{\left(\left(\mathbf{x}\_{i} - s\right)^{2} + z^{2}\right)\left(\mathbf{x}\_{i}\right)^{2}}{\left(\left(\mathbf{x}\_{i} - s\right)^{2} + z^{2}\right)^{2}}}{\left(\left(\mathbf{x}\_{i} - s\right)^{2} + z^{2}\right)^{2}} - 1$$

$$+2z\left[\frac{2\mathbf{x}\_i\mathbf{F}(4s^2+z^2)+12z^2\mathbf{s}}{\left(\mathbf{x}\_i^2+z^2\right)}-\frac{\left(\mathbf{x}\_i+\mathbf{s}\right)\mathbf{F}(4s^2+z^2)+6z^2\mathbf{s}}{\left(\mathbf{x}\_i+\mathbf{s}\right)^2+z^2}-\frac{\left(\mathbf{x}\_i-\mathbf{s}\right)\mathbf{F}(4s^2+z^2)+6z^2\mathbf{s}}{\left(\mathbf{x}\_i-\mathbf{s}\right)^2+z^2}\right]$$

Equation (22) can be solved for z using standard methods for solving nonlinear equations. Here, it is solved by an iterative method (Demidovich & Maron, 1973).

Substituting the computed depth (zc) as a fixed parameter in equation (16), we obtain

$$R(\mathbf{x}\_{i\prime}, \theta, \mathbf{s}) = \frac{z\_c(\mathbf{s}^2 + z\_c^2)R(\mathbf{0})}{2\mathbf{s}^2} P(\mathbf{x}\_{i\prime}, \theta, \mathbf{s}), \tag{23}$$

where

$$P(\mathbf{x}\_i, \theta, \mathbf{s}) = \left\{ \frac{2\mathbf{x}\_i \tan \theta + 2z\_c}{\mathbf{x}\_i^2 + z\_c^2} - \frac{(\mathbf{x}\_i + \mathbf{s})\tan \theta + z\_c}{\left(\mathbf{x}\_i + \mathbf{s}\right)^2 + z\_c^2} - \frac{(\mathbf{x}\_i - \mathbf{s})\tan \theta + z\_c}{\left(\mathbf{x}\_i - \mathbf{s}\right)^2 + z\_c^2} \right\}.$$

The unknown index parameter (θ) in equation (23) can be obtained by minimizing

$$\mathbb{V}\left(\theta\right) = \sum\_{i=1}^{N} \left[ L\left(\mathbf{x}\_{i}\right) - \frac{\mathbf{z}\_{c}\left(\mathbf{s}^{2} + \mathbf{z}\_{c}^{2}\right)R\left(\mathbf{0}\right)}{2\mathbf{s}^{2}} P(\mathbf{x}\_{i}, \theta, \mathbf{s})\right]^{2}. \tag{24}$$

Setting the derivative of ψ(θ) to zero with respect to θ leads to

$$\mathcal{A}(\theta) = \sum\_{i=1}^{N} \left[ L(\mathbf{x}\_i) - \frac{z\_c(\mathbf{s}^2 + z\_c^2)R(\mathbf{0})}{2\mathbf{s}^2} P(\mathbf{x}\_i, \theta, \mathbf{s}) \right] P^\*(\mathbf{x}\_i, \theta, \mathbf{s}) = \mathbf{0},\tag{25}$$

where P\*(xi, θ, s) is obtained by differentiating P(xi, θ, s) with respect to θ and is given by the following equation

$$P^\*(\mathbf{x}\_i, \theta, \mathbf{s}) = \text{Sec}^2 \theta \left[ \frac{2\mathbf{x}\_i}{\mathbf{x}\_i^2 + \mathbf{z}\_c^2} - \frac{(\mathbf{x}\_i + \mathbf{s})}{(\mathbf{x}\_i + \mathbf{s})^2 + \mathbf{z}\_c^2} - \frac{(\mathbf{x}\_i - \mathbf{s})}{(\mathbf{x}\_i - \mathbf{s})^2 + \mathbf{z}\_c^2} \right].$$

Equation (25) can be solved using a simple iterative method (Demodivch & Maron, 1973). Substituting the computed depth (zc) and index parameter (θc) in equation (14) as fixed parameters, we obtain

$$R(\mathbf{x}\_{i'}, \mathbf{s}) = KD(\mathbf{x}\_{i'}, \mathbf{s})\_{'} \tag{26}$$

where

30 New Achievements in Geoscience

where L(xi) denotes the moving average residual magnetic anomalies at xi. Setting the

\*

*f z Lx R W x zs W xzs*

*i i i*

( ) ( ) (0) ( , , ) ( , ) 0.

where W\*(xi, z, s) is evaluated by analytically differentiating W(xi, z, s) and is given by the

(( ) )(2( ) 12 ) 2 (( ) (4 ) 6 ) ( ,,) ( ) (( ) )

*x s z x s Fz zs z x s F s z z s W x zs s z*

*ii i*

2 (4 ) 12 ( ) (4 ) 6 ( ) (4 ) 6 <sup>2</sup> ( ) () ( )

*ii i*

Equation (22) can be solved for z using standard methods for solving nonlinear equations.

2 2 2 ( ) (0) ( , ,) ( , , ), <sup>2</sup> *c c i i zs z R Rx s Px s s*

2 tan 2 ( )tan ( )tan ( , ,) . () ()

( ) (0) () ( ) ( , ,) . <sup>2</sup> *<sup>N</sup> c c i i*

( ) (0) () ( ) ( , , ) ( , , ) 0, <sup>2</sup> *<sup>N</sup> c c <sup>i</sup> i i*

where P\*(xi, θ, s) is obtained by differentiating P(xi, θ, s) with respect to θ and is given by the

2 () () ( , ,) . () ()

2

*x xs xs P x s Sec*

*s*

*x z xs z xs z Px s*

The unknown index parameter (θ) in equation (23) can be obtained by minimizing

*ii i*

Substituting the computed depth (zc) as a fixed parameter in equation (16), we obtain

Here, it is solved by an iterative method (Demidovich & Maron, 1973).

1

 

*i*

Setting the derivative of ψ(θ) to zero with respect to θ leads to

 

1

\* 2

*i*

*i*

*i ii i*

*i ii*

2 2 22 2 2 22

( )(4 24 ) 2 (2 (4 ) 12 ) ( )

*x z x Fz zs z x F s z z s x z*

(22)

2 2 22 2

(( ) )(2( ) 12 ) 2 (( ) (4 ) 6 )

22 2 22 2 22 2 2 2 2 2 2 2

*xs z x s z x s Fz zs z x s F s z z s*

 2 2

(( )

*xF s z z s x sF s z z s x sF s z z s*

 

*x z xs z xs z*

2 2 2 2 2 2

*ic i c i c*

<sup>2</sup> 2 2

22 22 22

*ic i c i c*

*x z xs z xs z*

*ii i*

(25)

*x z xs z xs z*

2

2 2 \*

*zs z R L x Px s P x s*

*zs z R L x Px s s*

*i ci ci c*

*x s*

*i*

*i*

2 22 2 2 22 2 2

*z*

 

 

(24)

(23)

)

derivative of φ(z) to zero with respect to z leads to

following equation

*i*

where

\* 22

*z*

*i*

 

following equation

1

*i*

*N*

$$D(\mathbf{x}\_i, \mathbf{s}) = \frac{\mathbf{z}\_c}{2} \left\{ \frac{2\mathbf{x}\_i \sin \theta\_c + 2\mathbf{z}\_c \cos \theta\_c}{\mathbf{x}\_i^2 + \mathbf{z}\_c^2} - \frac{(\mathbf{x}\_i + \mathbf{s}) \sin \theta\_c + \mathbf{z}\_c \cos \theta\_c}{(\mathbf{x}\_i + \mathbf{s})^2 + \mathbf{z}\_c^2} - \frac{(\mathbf{x}\_i - \mathbf{s}) \sin \theta\_c + \mathbf{z}\_c \cos \theta\_c}{(\mathbf{x}\_i - \mathbf{s})^2 + \mathbf{z}\_c^2} \right\}.$$

Finally, applying the least-squares method to equation (26), the unknown amplitude coefficient (Kc) can be determined from

$$\mathbf{K}\_c = \frac{\sum\_{i=1}^{N} L(\mathbf{x}\_i) D(\mathbf{x}\_i, \mathbf{s})}{\sum\_{i=1}^{N} \left[ D(\mathbf{x}\_i, \mathbf{s}) \right]^2}. \tag{27}$$

Theoretically, one value of s is enough to determine the model parameters (z, θ, K) from equations (22), (25), and (27), respectively. In practice, more than one value of s is desirable because of the presence of noise, interference from neighbouring magnetic rocks, and complicated regional in data. Moreover, the advantage of this method over continuous modelling methods is that it requires neither magnetic susceptibility contrast nor depth information, and it can be applied if little or no factual information other than the magnetic data is available. The method may be automated to solve for a number of anomalies along a magnetic profile when the origin of each dike is known from other geophysical and/or geological data.

#### **2.2.1 Synthetic example**

The composite vertical magnetic anomaly in nanoteslas of Figure 6 consisting of the combined effect of an intermediate structure of interest (thin dike with K = 200 nT, z = 5 km, and θ = 40o), an interference neighbouring magnetic rocks (horizontal cylinder with K = 40 nT, z = 2 km, and θ = 30o), and a deep-seated structure (sphere with K = 80000 nT, z = 15 km, and θ = 10o) was computed by the following expression:

$$\Delta T(\mathbf{x}\_i) = 40 \frac{(4 - (\mathbf{x}\_i - 10)^2) \cos 30^\circ + 4(\mathbf{x}\_i - 10) \sin 30^\circ}{((\mathbf{x}\_i - 10)^2 + 2^2)^2}$$

$$\begin{aligned} & \quad \text{Horizontal cylinder (vertical component)}\\ & \quad + 200 \frac{\mathbf{x}\_i \sin 40^\circ + 5 \cos 40^\circ}{\mathbf{x}\_i^2 + 5^2} \\ & \quad \text{Thin dike (vertical component)}\\ & \quad + 80000 \frac{(450 - (\mathbf{x}\_i + 60)^2) \sin 10^\circ - 45(\mathbf{x}\_i + 60) \cos 10^\circ}{((\mathbf{x}\_i + 60)^2 + 15^2)^{\frac{5}{2}}} \\ & \quad \text{S sphere (vertical component)} \end{aligned} \tag{28}$$

New Techniques for Potential Field Data Interpretation 33

the buried dike. The result in this case is shown in Figure 7. The average depth obtained is 3.27 km compared with 4.99 km using the least-squares method. This result indicates that the present least-squares method is less sensitive to errors in the magnetic anomaly than

Table 4. Theoretical results obtained from the noisy composite magnetic anomaly. The depth, index parameter, and amplitude coefficient were computed from equations (22), (25), and (27), respectively. The actual dike parameters are: z = 5 km, θ = 40° and K = 200 nT.

Fig. 7. Euler deconvolution depth estimates from the data of the noisy composite vertical magnetic anomaly of Figure 6 using a thin dike model. The average depth obtained is 3.27 km.

anomaly profile of 24.64 m was digitised at an interval of 0.385 m.

sediments from the Parnaiba Basin, Brazil.

Figure 8 shows a total-field magnetic anomaly above a Mesozoic diabase dike (2 m thick) intruded into Palaeozoic sediments of the Parnaiba Basin, Brazil (Silva, 1989; Fig. 10). The total-field data were collected along profile perpendicular to the dike, with sampling interval of about 2 m. The depth to the outcropping dike (sensor height) is 1.9 m. This

Table 5. Computed values of depth, index parameter, and amplitude coefficient obtained from a total-field magnetic anomaly above a Mesozoic diabase dike intruded into Palaeozoic

Euler deconvolution method.

**2.2.2 Field example** 

Fig. 6. Noisy composite vertical magnetic anomaly consisting of the combined effect of an intermediate structure (thin dike with K = 200 nT, z = 5 km, and θ = 40°) [anomaly 1], a deep-seated structure (sphere with K = 80000 nT, z = 15 km, and θ = 10°) [anomaly 2], and an interference from neighbouring magnetic rocks (horizontal cylinder with K = 40 nT, z = 2 km, and θ = 30°) [anomaly 3].

In Figure 6, anomaly 1 is the anomaly due to the intermediate structure of our interest, anomaly 2 is the anomaly due to the interference from neighbouring magnetic rocks, and anomaly 3 is the anomaly due to the deep-seated structure. The composite magnetic anomaly ΔT(xi) is contaminated with 20% random error using the following equation:

$$
\Delta T\_{md}(\mathbf{x}\_i) = \Delta T(\mathbf{x}\_i) \left[ 1 + (RND(i) - 0.5) \, ^\circ 0.2 \right], \tag{29}
$$

where ΔTrnd(xi) is the contaminated anomaly value at xi, and RND(i) is a pseudo random number whose range is [0, 1]. The interval of the pseudo random number is an open interval, i.e. it does not include the extremes 0 and 1.

The noisy composite magnetic anomaly (ΔT) is subjected to a separation technique using the moving-average method. Nine successive moving-average graticule spacings were applied to each set of input data. Equations (22), (25), and (27) were applied to each of the nine moving-average residual profiles, yielding depth, index parameter, and amplitude coefficient solutions. The results and their average values are given in Table 4.

We verified numerically that the depth obtained is within 0.38%. The index parameter obtained is within 0.4%, whereas the error in the amplitude coefficient is within 4.7%. Good results are obtained by using the present algorithm for depth, index parameter, and amplitude coefficient, which are a primary concern in magnetic prospecting and other geophysical work. On the other hand, the Euler deconvolution method (Thompson, 1982) has been applied to the same noisy composite magnetic anomaly to determine the depth to

Fig. 6. Noisy composite vertical magnetic anomaly consisting of the combined effect of an intermediate structure (thin dike with K = 200 nT, z = 5 km, and θ = 40°) [anomaly 1], a deep-seated structure (sphere with K = 80000 nT, z = 15 km, and θ = 10°) [anomaly 2], and an interference from neighbouring magnetic rocks (horizontal cylinder with K = 40 nT, z = 2

In Figure 6, anomaly 1 is the anomaly due to the intermediate structure of our interest, anomaly 2 is the anomaly due to the interference from neighbouring magnetic rocks, and anomaly 3 is the anomaly due to the deep-seated structure. The composite magnetic anomaly ΔT(xi) is contaminated with 20% random error using the following equation:

where ΔTrnd(xi) is the contaminated anomaly value at xi, and RND(i) is a pseudo random number whose range is [0, 1]. The interval of the pseudo random number is an open

The noisy composite magnetic anomaly (ΔT) is subjected to a separation technique using the moving-average method. Nine successive moving-average graticule spacings were applied to each set of input data. Equations (22), (25), and (27) were applied to each of the nine moving-average residual profiles, yielding depth, index parameter, and amplitude

We verified numerically that the depth obtained is within 0.38%. The index parameter obtained is within 0.4%, whereas the error in the amplitude coefficient is within 4.7%. Good results are obtained by using the present algorithm for depth, index parameter, and amplitude coefficient, which are a primary concern in magnetic prospecting and other geophysical work. On the other hand, the Euler deconvolution method (Thompson, 1982) has been applied to the same noisy composite magnetic anomaly to determine the depth to

coefficient solutions. The results and their average values are given in Table 4.

*T x T x RND i rnd i* ( ) ( ) 1 ( ( ) 0.5) \* 0.2 , *<sup>i</sup>* (29)

km, and θ = 30°) [anomaly 3].

interval, i.e. it does not include the extremes 0 and 1.

the buried dike. The result in this case is shown in Figure 7. The average depth obtained is 3.27 km compared with 4.99 km using the least-squares method. This result indicates that the present least-squares method is less sensitive to errors in the magnetic anomaly than Euler deconvolution method.


Table 4. Theoretical results obtained from the noisy composite magnetic anomaly. The depth, index parameter, and amplitude coefficient were computed from equations (22), (25), and (27), respectively. The actual dike parameters are: z = 5 km, θ = 40° and K = 200 nT.

Fig. 7. Euler deconvolution depth estimates from the data of the noisy composite vertical magnetic anomaly of Figure 6 using a thin dike model. The average depth obtained is 3.27 km.

#### **2.2.2 Field example**

Figure 8 shows a total-field magnetic anomaly above a Mesozoic diabase dike (2 m thick) intruded into Palaeozoic sediments of the Parnaiba Basin, Brazil (Silva, 1989; Fig. 10). The total-field data were collected along profile perpendicular to the dike, with sampling interval of about 2 m. The depth to the outcropping dike (sensor height) is 1.9 m. This anomaly profile of 24.64 m was digitised at an interval of 0.385 m.


Table 5. Computed values of depth, index parameter, and amplitude coefficient obtained from a total-field magnetic anomaly above a Mesozoic diabase dike intruded into Palaeozoic sediments from the Parnaiba Basin, Brazil.

New Techniques for Potential Field Data Interpretation 35

**2.3 A window curves method to shape and depth solutions from third moving average** 

Following Abdelrahman & El-Araby (1993), the general gravity anomaly expression

2 2 ( ,,) , ( ) *<sup>i</sup> <sup>q</sup> i*

where q is the shape (shape factor), z is the depth, x is the position coordinate and A is amplitude coefficient related to the radius and density contrast of the buried structures. Examples of the shape factor for the semi-infinite vertical cylinder (3D), horizontal cylinder (2D), and sphere (3D) are 0.5, 1.0, and 1.5, respectively. Also, the shape factor for the finite vertical cylinder is approximately 1.0. The shape factor (q) approaches zero as the structure becomes a flat slab, and approaches 1.5 as the structure becomes a perfect sphere (point mass). The moving average (grid) method is an important yet simple technique for the separation of gravity anomalies into residual and regional components. The basic theory of the moving average methods is described by Griffin (1949), Agocs (1951), Abdelrahman &

Let us consider seven observation points (xi-3s, xi-2s, xi - s, xi, xi + s, xi +2s, and xi+3s) along the anomaly profile where s=1, 2,…, M spacing units and is called the window length. The first moving average regional gravity field Z1(xi, z, q, s) is defined as the average of g(xi-s, z,

1( , , ) (( ) ) (( ) ) . <sup>2</sup>

The first moving average residual gravity anomaly R1(xi, z, q, s) at the point xi is defined as

1( , , , ) 2( ) (( ) ) (( ) ) <sup>2</sup>

The second moving average residual gravity anomaly R2(xi, z, q, s) at the point xi is

( , , , ) 6( ) 4(( ) ) 4(( ) ) <sup>4</sup>

*i i*

2 2 2 2 2 2

6(( 2 ) ) 6(( 2 ) ) (( 3 ) ) (( 3 ) )

*xsz xsz xsz xsz*

2 2 2 2 2 2 2 2

*i ii i*

The third moving average residual gravity anomaly R3(xi, z, q, s) at any point xi is

( , , , ) 20( ) 15(( ) ) 15(( ) ) <sup>8</sup>

*i i*

 

*i i*

*i ii i*

For all shape factors (q), equation (34) gives the following value at xi=0

*<sup>A</sup> R x zqs x z x s z x s z*

*<sup>A</sup> R x zqs x z x s z x s z*

2 2 2 2

*q q*

*<sup>A</sup> Z x zs x s z x s z* (31)

2 2 2 2 2 2

2 2 2 2 2 2

(( 2 ) ) (( 2 ) )

*qqq*

*q q*

*q q*

 

*xsz xsz*

2 2 2 2

*qq q*

*q q*

(33)

(34)

*<sup>A</sup> R x zqs x z x s z x s z* (32)

*q qq*

q) and g(xi+s, z, q), which for the simple shapes mentioned above can be written as

*ii i*

*i ii i*

*x z*

*A*

(30)

produced by most geological structures is given by the following equation:

*gx zq*

**residual gravity anomalies** 

El-Araby (1996) and Abdelrahman et al. (2001).

(Abdelrahman & El-Araby, 1996).

(Abdelrahman et al. 2001)

3

2

Fig. 8. A total field magnetic anomaly above a Mesozoic diabase dike (2 m thick) intruded into Palaeozoic sediments from the Parnaiba Basin, Brazil (after Silva, 1989).

Fig. 9. Moving average residual magnetic anomalies above a Mesozoic diabase dike intruded into Palaeozoic sediments from the Parnaiba Basin, Brazil.

Adopting the same technique applied in the example above, we found the results given in Figure 9 and Table 5. The model parameters obtained are: z = 2.22 m, θ = 53o, and K = -215 nT. Silva (1989) applied a method to the same magnetic data employing the M-fitting technique and obtained a depth of about 3.5 m. The depth to the top is overestimated by both methods. This is reasonable because the upper part of the dike was weathered and the magnetite present was oxidized, losing most of its magnetic property (Silva, 1989). This field example shows also that there is a lack of "thinness". However, because many of the characteristics of wide dike anomalies are similar to those for thin dikes (Gay, 1963), our method can be applied not only to magnetic anomalies due to thin dikes, but also to those due to wide dikes to obtain reliable estimates of model parameters.

Fig. 8. A total field magnetic anomaly above a Mesozoic diabase dike (2 m thick) intruded

Fig. 9. Moving average residual magnetic anomalies above a Mesozoic diabase dike

Adopting the same technique applied in the example above, we found the results given in Figure 9 and Table 5. The model parameters obtained are: z = 2.22 m, θ = 53o, and K = -215 nT. Silva (1989) applied a method to the same magnetic data employing the M-fitting technique and obtained a depth of about 3.5 m. The depth to the top is overestimated by both methods. This is reasonable because the upper part of the dike was weathered and the magnetite present was oxidized, losing most of its magnetic property (Silva, 1989). This field example shows also that there is a lack of "thinness". However, because many of the characteristics of wide dike anomalies are similar to those for thin dikes (Gay, 1963), our method can be applied not only to magnetic anomalies due to thin dikes, but also to those

intruded into Palaeozoic sediments from the Parnaiba Basin, Brazil.

due to wide dikes to obtain reliable estimates of model parameters.

into Palaeozoic sediments from the Parnaiba Basin, Brazil (after Silva, 1989).

#### **2.3 A window curves method to shape and depth solutions from third moving average residual gravity anomalies**

Following Abdelrahman & El-Araby (1993), the general gravity anomaly expression produced by most geological structures is given by the following equation:

$$\log(\mathbf{x}\_i, z, \eta) = \frac{A}{(\mathbf{x}\_i^2 + z^2)^{\eta}},\tag{30}$$

where q is the shape (shape factor), z is the depth, x is the position coordinate and A is amplitude coefficient related to the radius and density contrast of the buried structures. Examples of the shape factor for the semi-infinite vertical cylinder (3D), horizontal cylinder (2D), and sphere (3D) are 0.5, 1.0, and 1.5, respectively. Also, the shape factor for the finite vertical cylinder is approximately 1.0. The shape factor (q) approaches zero as the structure becomes a flat slab, and approaches 1.5 as the structure becomes a perfect sphere (point mass). The moving average (grid) method is an important yet simple technique for the separation of gravity anomalies into residual and regional components. The basic theory of the moving average methods is described by Griffin (1949), Agocs (1951), Abdelrahman & El-Araby (1996) and Abdelrahman et al. (2001).

Let us consider seven observation points (xi-3s, xi-2s, xi - s, xi, xi + s, xi +2s, and xi+3s) along the anomaly profile where s=1, 2,…, M spacing units and is called the window length. The first moving average regional gravity field Z1(xi, z, q, s) is defined as the average of g(xi-s, z, q) and g(xi+s, z, q), which for the simple shapes mentioned above can be written as

$$Z\_1(\mathbf{x}\_i, z, \mathbf{s}) = \frac{A}{2} \left\{ (\left(\mathbf{x}\_i - \mathbf{s}\right)^2 + z^2)^{-q} + \left(\left(\mathbf{x}\_i + \mathbf{s}\right)^2 + z^2\right)^{-q} \right\}.\tag{31}$$

The first moving average residual gravity anomaly R1(xi, z, q, s) at the point xi is defined as (Abdelrahman & El-Araby, 1996).

$$R\_1(\mathbf{x}\_i, z, q, s) = \frac{A}{2} \left[ \mathbf{2} (\mathbf{x}\_i^2 + z^2)^{-q} - \left( (\{\mathbf{x}\_i - s\}^2 + z^2)^{-q} + (\{\mathbf{x}\_i + s\}^2 + z^2)^{-q} \right) \right].\tag{32}$$

The second moving average residual gravity anomaly R2(xi, z, q, s) at the point xi is (Abdelrahman et al. 2001)

$$\begin{split} R\_2(\mathbf{x}\_i, z, q\_i, \mathbf{s}) &= \frac{A}{4} \Big[ 6\left(\mathbf{x}\_i^2 + z^2\right)^{-q} - 4\left(\left(\mathbf{x}\_i - \mathbf{s}\right)^2 + z^2\right)^{-q} - 4\left(\left(\mathbf{x}\_i + \mathbf{s}\right)^2 + z^2\right)^{-q} \\ &\quad + \left(\left(\mathbf{x}\_i - 2\mathbf{s}\right)^2 + z^2\right)^{-q} + \left(\left(\mathbf{x}\_i + 2\mathbf{s}\right)^2 + z^2\right)^{-q} \Big]. \end{split} \tag{33}$$

The third moving average residual gravity anomaly R3(xi, z, q, s) at any point xi is

$$R\_3(\mathbf{x}\_i, z, q, s) = \frac{A}{8} \Big[ 20(\mathbf{x}\_i^2 + z^2)^{-q} - 15(\left(\mathbf{x}\_i + s\right)^2 + z^2)^{-q} - 15(\left(\mathbf{x}\_i - s\right)^2 + z^2)^{-q}$$

$$+ \mathsf{6}(\left(\mathbf{x}\_i + 2s\right)^2 + z^2)^{-q} + \mathsf{6}(\left(\mathbf{x}\_i - 2s\right)^2 + z^2)^{-q} \tag{34}$$

$$- \left(\left(\mathbf{x}\_i + 3s\right)^2 + z^2\right)^{-q} - \left(\left(\mathbf{x}\_i - 3s\right)^2 + z^2\right)^{-q} \Big].$$

For all shape factors (q), equation (34) gives the following value at xi=0

New Techniques for Potential Field Data Interpretation 37

each window length, a depth value was determined iteratively for all shape values, and a

We have also applied the first moving average method (Abdelrahman & El-Araby 1996) and the second moving average method (Abdelrahman et al. 2001) to the same input data generated from equation (38). The results in the case of using first, second, and third moving

Fig. 10. Composite gravity anomaly (Δg) of a buried horizontal cylinder and a dipping fault

Fig. 11. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the first moving average method. Estimates of q and z are,

window curve was plotted, illustrating the relation between the depth and shape.

average techniques are shown in Figures 11, 12, and 13, respectively.

as obtained from equation (38).

respectively, 0.26, and 2.1 km.

$$A = \frac{4R\_3(0)}{\left[10z^{-2q} - 15(s^2 + z^2)^{-q} + 6(4s^2 + z^2)^{-q} - (9s^2 + z^2)^{-q}\right]} \tag{35}$$

Using equation (34) and equation (35), we obtain the following normalized equation at xi=s

$$\frac{R\_3(\mathbf{s})}{R\_3(0)} = \frac{\left[\frac{26(\mathbf{s}^2 + \mathbf{z}^2)^{-q} - 16(4\mathbf{s}^2 + \mathbf{z}^2)^{-q} + 6(9\mathbf{s}^2 + \mathbf{z}^2)^{-q} - 15\mathbf{z}^{-2q} - (16\mathbf{s}^2 + \mathbf{z}^2)^{-q}\right]}{\left[\frac{20}{2}\mathbf{z}^{-2q} - 30(\mathbf{s}^2 + \mathbf{z}^2)^{-q} + 12(4\mathbf{s}^2 + \mathbf{z}^2)^{-q} - 2(9\mathbf{s}^2 + \mathbf{z}^2)^{-q}\right]} \text{.} \tag{36}$$

Equation (36) is independent on the amplitude coefficient (A). Let F= R3(s)/R3(0) then from equation (36) we obtain

$$z = \left[\frac{-15}{F \ast P(s, z, q) + \mathcal{W}(s, z, q)}\right]^{\frac{1}{2\gamma}},\tag{37}$$

Where

$$P(s,z,q) = \left[20z^{-2q} - 30(s^2 + z^2)^{-q} + 12(4s^2 + z^2)^{-q} - 2(9s^2 + z^2)^{-q}\right]\_{s'}$$

and

$$\mathcal{W}(\mathbf{s}, z, \mathbf{q}) = \left[ 16(\mathbf{4}\mathbf{s}^2 + z^2)^{-q} - 6(\mathbf{9}\mathbf{s}^2 + z^2)^{-q} + (\mathbf{16}\mathbf{s}^2 + z^2)^{-q} - 2\mathbf{6}(\mathbf{s}^2 + z^2)^{-q} \right].$$

Equation (37) can be used not only to determine the depth, but also to simultaneously estimate the shape of the buried structure as will be illustrated by theoretical and field examples.

#### **2.3.1 Synthetic example**

The composite gravity anomaly in mGals of Figure 10 consisting of the combined effect of a shallow structure (horizontal cylinder with 3 km depth) and a deep-seated structure (dipping fault with depth to upper faulted slab = 15 km, depth to lower faulted slab = 20 km, and dip angle of fault plane = 50 degrees) was computed by the above expression:

$$\Delta \mathbf{g}(\mathbf{x}\_i) = \frac{200}{\left(\mathbf{x}\_i^2 + \mathbf{3}^2\right)} \quad + \quad 100(\tan^{-1}(\frac{\mathbf{x}\_i - \mathbf{5}}{15} + \cot(50^\circ)) - \tan^{-1}(\frac{\mathbf{x}\_i - \mathbf{5}}{20} + \cot(50^\circ))), \tag{38}$$

(Horizontal cylinder + dipping fault)

The composite gravity field (Δg) was subjected to our third moving average technique. The third moving average residual value at point xi is computed from the input gravity data g(xi) using the following equation:

$$R\_3(\mathbf{x}\_i) = \frac{20\mathbf{g}(\mathbf{x}\_i) - 15\mathbf{g}(\mathbf{x}\_i + \mathbf{s}) - 15\mathbf{g}(\mathbf{x}\_i - \mathbf{s}) + 6\mathbf{g}(\mathbf{x}\_i - 2\mathbf{s}) + 6\mathbf{g}(\mathbf{x}\_i + 2\mathbf{s}) - \mathbf{g}(\mathbf{x}\_i + 3\mathbf{s}) - \mathbf{g}(\mathbf{x}\_i - 3\mathbf{s})}{8}. \tag{39}$$

Five successive third moving average windows (s= 2, 3, 4, 5, and 6 km) were applied to input data. Each of the resulting moving average profiles was analyzed by equation (37). For

3 2 22 22 22 4 (0) 10 15( ) 6(4 ) (9 ) *q q qq*

(35)

(36)

(37)

  (38)

*z sz sz sz*

22 22 22 2 2 2

*q q qq q*

1

*q*

2 22 22 22

<sup>2</sup> <sup>15</sup> , \* (, , ) (, , )

 

*qq qq*

Using equation (34) and equation (35), we obtain the following normalized equation at xi=s

26( ) 16(4 ) 6(9 ) 15 (16 ) ( )

Equation (36) is independent on the amplitude coefficient (A). Let F= R3(s)/R3(0) then from

*F Pszq Wszq* 

<sup>2</sup> 22 22 22 ( , , ) 20 30( ) 12(4 ) 2(9 ) , *qq qq Pszq z s z s z s z*

22 22 22 22 ( , , ) 16(4 ) 6(9 ) (16 ) 26( ) . *qqqq Wszq s z s z s z s z*

Equation (37) can be used not only to determine the depth, but also to simultaneously estimate the shape of the buried structure as will be illustrated by theoretical and field

The composite gravity anomaly in mGals of Figure 10 consisting of the combined effect of a shallow structure (horizontal cylinder with 3 km depth) and a deep-seated structure (dipping fault with depth to upper faulted slab = 15 km, depth to lower faulted slab = 20 km, and dip angle of fault plane = 50 degrees) was computed by the above expression:

<sup>200</sup> 5 5 ( ) 100(tan ( cot(50 )) tan ( cot(50 ))), ( 3) <sup>15</sup> <sup>20</sup>

The composite gravity field (Δg) was subjected to our third moving average technique. The third moving average residual value at point xi is computed from the input gravity data

20 ( ) 15 ( ) 15 ( ) 6 ( 2 ) 6 ( 2 ) ( 3 ) ( 3 ) ( ) . <sup>8</sup> *ii i i i i i*

Five successive third moving average windows (s= 2, 3, 4, 5, and 6 km) were applied to input data. Each of the resulting moving average profiles was analyzed by equation (37). For

*<sup>g</sup> <sup>x</sup> <sup>g</sup> x s <sup>g</sup> x s <sup>g</sup> x s <sup>g</sup> x s <sup>g</sup> x s <sup>g</sup> x s R x* (39)

1 1

*i i*

*x x*

 

(0) 20 30( ) 12(4 ) 2(9 )

*R z sz sz sz*

*R s sz sz sz z sz*

*<sup>R</sup> <sup>A</sup>*

*z*

3

3

equation (36) we obtain

Where

and

examples.

3

*i*

**2.3.1 Synthetic example** 

*i*

*g x*

2 2

*i*

*x*

(Horizontal cylinder + dipping fault)

g(xi) using the following equation:

each window length, a depth value was determined iteratively for all shape values, and a window curve was plotted, illustrating the relation between the depth and shape.

We have also applied the first moving average method (Abdelrahman & El-Araby 1996) and the second moving average method (Abdelrahman et al. 2001) to the same input data generated from equation (38). The results in the case of using first, second, and third moving average techniques are shown in Figures 11, 12, and 13, respectively.

Fig. 10. Composite gravity anomaly (Δg) of a buried horizontal cylinder and a dipping fault as obtained from equation (38).

Fig. 11. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the first moving average method. Estimates of q and z are, respectively, 0.26, and 2.1 km.

New Techniques for Potential Field Data Interpretation 39

whereas the model parameters obtained from the third moving average method are in excellent agreement with the parameters of the shallow structure given in model [equation (38)]. The conclusion is inescapable that the third moving average technique gives the best result. Moreover, random errors were added to the composite gravity anomaly Δg(xi) to produce noisy anomaly. In this case, we choose a white random noise with amplitude being 5% of the maximum amplitude variation along the entire profile, so that the noise will be equally along the profile. The noisy anomaly is subjected to a separation technique using only the third moving average method. Five successive third moving average windows were applied to the noisy input data. Adapting the same

Fig. 14. Family of window curves of z as a function of q for s= 2, 3, 4, 5, 6, and 7 km as obtained from gravity anomaly (Δg) after adding 5% random errors to the data using the present third

Figure 14 shows the window curves intersect at approximately q= 1.09 and z= 2.98 km. The results (Fig. 14) are generally in very good agreement with the shape factor and depth parameters shown in [equation (38)]. This indicates that our method will produce reliable

A Bouguer gravity anomaly profile along AÁ of the gravity map of the Humble salt dome, near Houston (Nettleton, 1962, Fig. 22) is shown in Figure 15. The gravity profile has been digitized at an interval of 0.26 km. The Bouguer gravity anomalies thus obtained have been subjected to a separation technique using the third moving average method. Filters were applied in four successive windows (s= 2.60, 2.86, 3.12, and 3.38 km). In this way, four third moving average residual anomaly profiles were obtained (Fig. 16). The same procedure described for the synthetic examples was used to estimate the shape and the depth of the salt dome. The results are plotted in Figure 17. Figure 17 shows that the window curves intersect at a point where the depth is about 4.95 km and at q=1.41. This

moving average method. Estimates of q and z are, respectively, 1.09, and 2.98 km.

shape and depth estimates when applied to noisy gravity data.

**2.3.2 Field example** 

interpretation procedure, the results are shown in Figure 14.

Fig. 12. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the second moving average method. Estimates of q and z are, respectively, 0.75, and 2.7 km.

Fig. 13. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the present third moving average method. Estimates of q and z are, respectively, 1.03, and 3.05 km.

In the case of using first moving average technique, the window curves intersect at a narrow region where 0.55>q>0.05 and 2.5 km>z> 1.8 km (Fig. 11). The central point of this region occurs at the location q= 0.26 and z= 2.1 km. On the other hand, in the case of using second moving average technique, the curves intersect each other nearly at a point where q= 0.75 and z= 2.7 km (Fig. 12). Finally, in the case of using the third moving average method (present method), the curves intersect each other at the correct locations q= 1.03 and z= 3.05 km (Fig. 13). From these results (Figs. 11-13), it is evident that the first and the second moving average methods give erroneous results for the shape and depth solutions

Fig. 12. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the second moving average method. Estimates of q and z

Fig. 13. Family of window curves of z as a function of q for s= 2, 3, 4, 5, and 6 km as obtained from gravity anomaly (Δg) using the present third moving average method. Estimates of q

In the case of using first moving average technique, the window curves intersect at a narrow region where 0.55>q>0.05 and 2.5 km>z> 1.8 km (Fig. 11). The central point of this region occurs at the location q= 0.26 and z= 2.1 km. On the other hand, in the case of using second moving average technique, the curves intersect each other nearly at a point where q= 0.75 and z= 2.7 km (Fig. 12). Finally, in the case of using the third moving average method (present method), the curves intersect each other at the correct locations q= 1.03 and z= 3.05 km (Fig. 13). From these results (Figs. 11-13), it is evident that the first and the second moving average methods give erroneous results for the shape and depth solutions

are, respectively, 0.75, and 2.7 km.

and z are, respectively, 1.03, and 3.05 km.

whereas the model parameters obtained from the third moving average method are in excellent agreement with the parameters of the shallow structure given in model [equation (38)]. The conclusion is inescapable that the third moving average technique gives the best result. Moreover, random errors were added to the composite gravity anomaly Δg(xi) to produce noisy anomaly. In this case, we choose a white random noise with amplitude being 5% of the maximum amplitude variation along the entire profile, so that the noise will be equally along the profile. The noisy anomaly is subjected to a separation technique using only the third moving average method. Five successive third moving average windows were applied to the noisy input data. Adapting the same interpretation procedure, the results are shown in Figure 14.

Fig. 14. Family of window curves of z as a function of q for s= 2, 3, 4, 5, 6, and 7 km as obtained from gravity anomaly (Δg) after adding 5% random errors to the data using the present third moving average method. Estimates of q and z are, respectively, 1.09, and 2.98 km.

Figure 14 shows the window curves intersect at approximately q= 1.09 and z= 2.98 km. The results (Fig. 14) are generally in very good agreement with the shape factor and depth parameters shown in [equation (38)]. This indicates that our method will produce reliable shape and depth estimates when applied to noisy gravity data.

#### **2.3.2 Field example**

A Bouguer gravity anomaly profile along AÁ of the gravity map of the Humble salt dome, near Houston (Nettleton, 1962, Fig. 22) is shown in Figure 15. The gravity profile has been digitized at an interval of 0.26 km. The Bouguer gravity anomalies thus obtained have been subjected to a separation technique using the third moving average method. Filters were applied in four successive windows (s= 2.60, 2.86, 3.12, and 3.38 km). In this way, four third moving average residual anomaly profiles were obtained (Fig. 16). The same procedure described for the synthetic examples was used to estimate the shape and the depth of the salt dome. The results are plotted in Figure 17. Figure 17 shows that the window curves intersect at a point where the depth is about 4.95 km and at q=1.41. This

New Techniques for Potential Field Data Interpretation 41

Fig. 17. Family of window curves of z as a function of q for s= 2.60, 2. 86, 3.12 and 3.38 km as

Three different numerical algorithms to depth and other associated model parameters determination from magnetic and gravity data have been presented. In the first method, a least-squares minimization approach has been developed to determine the depth to a buried structure from magnetic data. The problem of determining the depth of a buried structure from the magnetic anomaly has been transformed into the problem of finding a solution of a nonlinear equation. The method presented is very simple to execute. The advantages of the present method over previous techniques, which use only a few points, distances, standardized curves, and nomograms are: (1) all observed values cal be used, (2) the method is automatic, and (3) the method is not sensitive to errors in the magnetic anomaly. Finally, the advantage of the proposed method over previous least-squares techniques is that any initial guess for the depth parameter works well. In the second method, the problem of determining the depth, index parameter, and amplitude coefficient of the buried thin dike from a magnetic anomaly profile has been transformed into the problem of solving two individual nonlinear equations and a linear equation, respectively. A scheme for interpreting the magnetic data to obtain the model parameters based on the least-squares method is developed. The method provides two advantages over previous least-squares techniques: (1) each model parameter is computed from all observed data, and (2) the method is less sensitive to errors in the magnetic anomaly. The method involves using a thin dike model convolved with the same moving average filter as applied to the observed data. As a result, the method can be applied not only to residuals, but also to measured magnetic data having noise, interference from neighboring magnetic rocks, and complicated regional. Synthetic and field studies demonstrate the efficiency of the present technique. Finally, the problem of determining the shape and depth from gravity data of long or short profile length can be solved using the window curves methods for simple anomalies. The window curves methods are very simple to execute and work well even when the data contains

obtained from the Humble gravity anomaly profile

**3. Conclusion** 

suggests that the shape of the salt dome resembles a sphere or, practically, a threedimensional source with a hemi- spherical roof and root. This result is generally in good agreement with the dome form estimated from drilling top and contact (Nettleton, 1976; Fig. 8-16).

Fig. 15. Observed gravity profile of the Humble salt dome, near Houston, TX, USA (after Nettleton 1962)

Fig. 16. Third moving average residual gravity anomalies on line AÁ of the Humble salt dome, for s= 2.60, 2. 86, 3.12 and 3.38 km.

Fig. 17. Family of window curves of z as a function of q for s= 2.60, 2. 86, 3.12 and 3.38 km as obtained from the Humble gravity anomaly profile

#### **3. Conclusion**

40 New Achievements in Geoscience

suggests that the shape of the salt dome resembles a sphere or, practically, a threedimensional source with a hemi- spherical roof and root. This result is generally in good agreement with the dome form estimated from drilling top and contact (Nettleton, 1976;

Fig. 15. Observed gravity profile of the Humble salt dome, near Houston, TX, USA (after

Fig. 16. Third moving average residual gravity anomalies on line AÁ of the Humble salt

dome, for s= 2.60, 2. 86, 3.12 and 3.38 km.

Fig. 8-16).

Nettleton 1962)

Three different numerical algorithms to depth and other associated model parameters determination from magnetic and gravity data have been presented. In the first method, a least-squares minimization approach has been developed to determine the depth to a buried structure from magnetic data. The problem of determining the depth of a buried structure from the magnetic anomaly has been transformed into the problem of finding a solution of a nonlinear equation. The method presented is very simple to execute. The advantages of the present method over previous techniques, which use only a few points, distances, standardized curves, and nomograms are: (1) all observed values cal be used, (2) the method is automatic, and (3) the method is not sensitive to errors in the magnetic anomaly. Finally, the advantage of the proposed method over previous least-squares techniques is that any initial guess for the depth parameter works well. In the second method, the problem of determining the depth, index parameter, and amplitude coefficient of the buried thin dike from a magnetic anomaly profile has been transformed into the problem of solving two individual nonlinear equations and a linear equation, respectively. A scheme for interpreting the magnetic data to obtain the model parameters based on the least-squares method is developed. The method provides two advantages over previous least-squares techniques: (1) each model parameter is computed from all observed data, and (2) the method is less sensitive to errors in the magnetic anomaly. The method involves using a thin dike model convolved with the same moving average filter as applied to the observed data. As a result, the method can be applied not only to residuals, but also to measured magnetic data having noise, interference from neighboring magnetic rocks, and complicated regional. Synthetic and field studies demonstrate the efficiency of the present technique. Finally, the problem of determining the shape and depth from gravity data of long or short profile length can be solved using the window curves methods for simple anomalies. The window curves methods are very simple to execute and work well even when the data contains

New Techniques for Potential Field Data Interpretation 43

Chai, Y., & Hinze, W. J. (1988). Gravity inversion of an interface above which the density

Demidovich, B., & Maron, I. A. (1973). *Computational mathematics,* MIR Publishers, ISBN

Essa, K. S. (2007), A simple formula for shape and depth determination from residual gravity anomalies. *Acta Geophysica,* Vol. 55, No. 2, pp. 182-190, ISSN 1895-6572. Essa, K. S. (2011), A new algorithm for gravity or self-potential data interpretation. *Journal of Geophysics and Engineering,* Vol. 8, No. 3, pp. 434-446, ISSN 1742-2132. Gay, P. (1963). Standard curves for interpretation of magnetic anomalies over long tabular

Griffin, W. R. (1949). Residual gravity in theory and practice. *Geophysics,* Vol. 14, No. 1, pp.

Keating, P., & Pilkington, M. (2004). Euler deconvolution of the analytic signal and its

Martin-Atienza, B., & Garcia-Abdeslem, J. (1999). 2-D gravity modeling with analytically

Martins, C. M., Barbosa, V. C. F., & Silva, J. B. C. (2010 ). Simultaneous 3D depth-to-

McGrath, P. H., & Hood, P. J. (1973). An automatic least-squares multimodel method for magnetic interpretation. *Geophysics,* Vol. 38, No. 2, pp. 349-358, ISSN 0016-8033. Mohan, N. L., Sundararajan, N. & Seshagiri Rao, S. V. (1982). Interpretation of some two-

Nabighian, M. N. (1984). Toward a three-dimensional automatic interpretation of potential

Nettleton, L. L. (1962). Gravity and magnetic for geologists and seismologists. AAPG, Vol.

Nettleton, L. L. (1976). *Gravity and magnetics in oil prospecting,* Mc-Graw Hill Book Co., ISBN

Nunes, T. M., Barbosa, V. C. F., & Silva, J. B. C. (2008). Magnetic basement depth inversion

Prakasa Rao, T. K. S. & Subrahmanyam, M. (1988). Characteristic curves for the inversion of

Prakasa Rao, T. K. S., Subrahmanyam, M., & Srikrishna Murthy, A. (1986). Nomogram for

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. (1986). *Numerical recipes,* 

*Geophysics,* Vol. 51, No. 11, pp. 2156-2159, ISSN 0016-8033.

few points. *Geophysics,* Vol. 75, No. 3, pp. I21–I28, ISSN 0016-8033.

application to magnetic interpretation. Geophysical Prospecting, Vol. 52, No. 3, pp.

defined geometry and quadratic polynomial density functions. *Geophysics,* Vol. 64,

basement and density-contrast estimates using gravity data and depth control at

dimensional magnetic bodies using Hilbert transform. *Geophysics,* Vol. 47, No. 3,

field data via generalized Hilbert transforms: Fundamental relations. *Geophysics,*

in the space domain. *Pure Applied Geophysics,* Vol. 165, No. 9-10, pp. 891–1911, ISSN

magnetic anomalies of spherical ore bodies. *Pure and Applied Geophysics,* Vol. 126,

the direct interpretation of magnetic anomalies due to long horizontal cylinders.

*the art of scientific computing,* Cambridge University Press, ISBN 0-521-43064-X,

bodies. *Geophysics,* Vol. 28, No. 2, pp. 161-200, ISSN 0016-8033.

ISSN 0016-8033.

B000R05KIU, Moscow.

39-56, ISSN 0016-8033.

165–182, ISSN 1365-2478 .

pp. 376-387, ISSN 0016-8033.

9780070463035, New York.

No. 1, pp. 69-83, ISSN 0033-4553.

0033-4553.

Cambridge.

Vol. 49, No. 6, pp. 780-786, ISSN 0016-8033.

46, No. , pp. 1815-1838, ISSN 0149-1423.

No. 6, pp. 1730–1734, ISSN 0016-8033.

contrast varies exponentially with depth. *Geophysics,* Vol. 53, No. 6, pp. 837–845,

noise. The methods involve using simple models convolved with the same moving average filter as applied to the observed gravity data. As a result, the methods can be applied not only to residuals, but also to measured gravity data. Given their relative strength, the present techniques complement the existing methods of shape and depth determination from moving average residual gravity anomaly profiles and overcome some of their shortcomings. The present iterative methods converge to a depth solution for any shape factor value and when applying any window length. The methods are developed to derive shape–related parameter and depth and remove the regional trend from gravity data. These two parameters might be used to gain geological insight concerning the subsurface. However, the advantages of the third moving average method over the first and the second moving average techniques are: (1) the method removes adequately the regional field due to deep-seated structure from gravity data and (2) the method can be applied to large gridded data. Theoretical and field examples have illustrated the validity of the algorithms presented.

## **4. Acknowledgment**

The author thanks Prof. Aleksander Lazinica, CEO, Prof. Hwee-San Lim, the Book editor and Mrs. Ana Skalamera, Publishing Process Manager for their excellent suggestions and collaborations. He would like to thank Prof. Sharf El Din Mahmoud, Chairman and Prof. El-Sayed Abdelrahman, Geophysics Department, Faulty of Science, Cairo University for their continuous support and encouragement.

#### **5. References**


noise. The methods involve using simple models convolved with the same moving average filter as applied to the observed gravity data. As a result, the methods can be applied not only to residuals, but also to measured gravity data. Given their relative strength, the present techniques complement the existing methods of shape and depth determination from moving average residual gravity anomaly profiles and overcome some of their shortcomings. The present iterative methods converge to a depth solution for any shape factor value and when applying any window length. The methods are developed to derive shape–related parameter and depth and remove the regional trend from gravity data. These two parameters might be used to gain geological insight concerning the subsurface. However, the advantages of the third moving average method over the first and the second moving average techniques are: (1) the method removes adequately the regional field due to deep-seated structure from gravity data and (2) the method can be applied to large gridded data. Theoretical and field examples have illustrated the validity of the algorithms

The author thanks Prof. Aleksander Lazinica, CEO, Prof. Hwee-San Lim, the Book editor and Mrs. Ana Skalamera, Publishing Process Manager for their excellent suggestions and collaborations. He would like to thank Prof. Sharf El Din Mahmoud, Chairman and Prof. El-Sayed Abdelrahman, Geophysics Department, Faulty of Science, Cairo University for their

Abdelrahman, E. M., & El-Araby, H. M. (1993). Shape and depth solutions from gravity data

Abdelrahman, E. M., & El-Araby, T. M. (1996). Shape and depth solutions from moving

Abdelrahman, E. M., El-Araby, T. M., El-Araby, H. M., & Abo-Ezz, E. R. (2001). A new

Abdelrahman, E. M., & Essa, K. S. (2005). Magnetic interpretation using a least-squares,

Agocs, W. B. (1951). Least squares residual anomaly determination. *Geophysics,* Vol. 16, No.

Atchuta Rao, D. A., & Ram Babu, H. V. (1980). Properties of the Relation Figures between

Atchuta Rao, D., Ram Babu, H. V., & Sankar Narayan, P. V. (1980). Relationship of magnetic

*Geophysics,* Vol. 45, No. 1, pp. 32-36, ISSN 0016-8033.

using correlation factors between successive least-squares residuals. *Geophysics,*

average residual gravity anomalies. *Journal of Applied Geophysics,* Vol. 36, No. 2-3,

method for shape and depth determinations from gravity data. *Geophysics,* Vol. 66,

depth-shape curves method. *Geophysics,* Vol. 70, No. 3, pp. L23–L40, ISSN 0016-

the Total, Vertical, and Horizontal Field Magnetic Anomalies over a Long Horizontal Cylinder Ore Body. *Current Science,* Vol. 49, No. 15, pp. 584–585, ISSN

anomalies due to subsurface features and the interpretation of sloping contacts.

presented.

**4. Acknowledgment** 

**5. References** 

8033.

001-3891.

continuous support and encouragement.

pp. 89-95, ISSN 0926-9851.

No. 6, pp. 1774-1780, ISSN 0016-8033.

4, pp. 686-696, ISSN 0016-8033.

Vol. 58, No. 12, pp. 1785-1791, ISSN 0016-8033.


**3** 

*Poland* 

Jadwiga A. Jarzyna et al.\*

*AGH University of Science and Technology,* 

**Geophysics in Near-Surface Investigations** 

Environmental and engineering geophysics are new branches of geophysics that focus on data collection and an analysis of the present day environment. Noninvasive, nondestructive, and inexpensive geophysical technologies are now capable of providing detailed 3D and 4D representations of the subsurface. These advances in geophysical survey techniques have not only enhanced traditional applications, such as prospecting geophysics, but also made way for development of new methodologies in near surface investigations. Environmental and engineering geophysics includes mapping of shallow geological formations, prospecting for resources located in the near surface, monitoring aquifers, mapping anthropogenic effects on the environment, and surveying for civil engineering and archaeological projects. Environmental and engineering geophysics measures parameters such as salinity of subsurface fluids, levels of radioactivity, and other soil properties affected

The increasing demand for sustainable development as well as community efforts to preserve and conserve the environment have fostered the development of engineering technologies that are more environmentally friendly than traditional geophysical tools. Newer techniques are often used to monitor and protect the environment, as well as evaluate geotechnical risks of various natural and human - made structures. The increasing number of these modern geophysical applications, and of professionals specializing in them,

**2. Petrophysical background for geophysical investigations (Jadwiga A.** 

*AGH University of Science and Technology, Faculty of Geology Geophysics and Environmental Protection,* 

Jerzy Dec, Jerzy Karczewski, Sławomir Porzucek, Sylwia Tomecka-Suchoń,

Petrophysics is a subdiscipline of geophysics that addresses the physical properties and other parameters of rock and other subsurface material. Applied geophysics combines the theoretical foundation of a given geophysical method with empirical understanding of the materials involved. In the case of petrophysics, theoretical parameters of rock and other subsurface materials (i.e. their physical properties) are used to interpret empirical

**1. Introduction** 

by industrial activity.

**Jarzyna)** 

*Krakow, Poland* 

Anna Wojas and Jerzy Ziętek

 \*

bear witness to their emerging significance.

*Faculty of Geology Geophysics and Environmental Protection, Krakow,* 

