3D Modeling and Inversion of Gravity Data in Exploration Scale 3D Modeling and Inversion of Gravity Data in Exploration Scale

The gravity method has been widely used for detecting the subsurface density anomaly and geological structures. The interpretation result based on gravity data can be used for mineral/oil exploration and regional geological study. The effective and successful appli- cation of gravity methods depends on the fast forward modeling and stable inversion tools to image the subsurface density structures. In this chapter, we will review the applications and developments of gravity method. We starts form the basic theory for gravity field and the scalar gravity potential and introduce the closed form of the solution for the gravity field caused by a density anomaly. Different gravity data forward modeling and inversion techniques will be introduced in this chapter with their application in petroleum reconnaissance. Several examples will be presented in this chapter to illustrate the application of different gravity modeling and inversion techniques.


Introduction
Gravity method is an important geophysical tool for subsurface mapping of the density distributions. The development of high-accuracy instrumentation and stable observation platform since 1980s has made the large scale airborne gravity survey become possible. The interpretation result based on gravity data can be used for mineral and oil exploration. For example, the gravity data has been used to estimate the depth of the sedimentary basin and detect sub-sea basalt which is usually associated with the oil reservoir.
The gravity modeling has its origin back to Newton's law of universal gravitation developed in 1687, which is equivalent to solving a Poisson's equation for scalar gravity potential. One can use different approaches to solve this Poisson's equation, e.g., finite difference or finite element method. However, the most common approach is to solve it in an integral way using Green's function. In the first section of this chapter, we will review the Poisson's equation and its integral solution for a general three-dimensional case. The numerical approaches to evaluate the integrals in the solution of the Poisson's equation will be discussed. For more realistic problems, we also provide an approach of gravity modeling to take into account of complex geometry such as topography.
Conventional gravity inversion is based on discretizing the earth model into a set of rectangular prisms, and considering each prism has a constant density. Given observed gravity and/or gravity gradiometry fields, one is interested to recover the density within each cell. Similar as other geophysical inverse problems, this volumetric inversion of gravity data is an ill-posed problem which is strongly affected by the non-uniqueness and instability. To solve this illposed a problem, we will first introduce an objective functional which contains two parts, the misfit functional for least-square fitting of the observed data, and the stabilizer for regularizing the subsurface density distributions. The solutions for solving the minimization problem of this objective functional will then be briefly discussed. Generally, the volumetric gravity inversion produces diffusive images of the subsurface density distribution with very limited depth resolution (strictly speaking, the gravity method, with observation above the earth's surface, does not have any depth resolution for the subsurface density distribution). One may improve the gravity inversion results by adding more constraints, e.g., logarithmic barrier to constrain upper and lower limits of the density, focusing stabilizer to produce sharper boundaries between the density contrasts, and etc., in the inversion.
A more advanced approach of gravity inversion is to invert for the shape of the anomaly instead of the three-dimensional volumetric density distribution. In certain applications, such as sedimentary basin analysis, it is more interesting to estimate the location of the sedimentbasement interface and the shape by assuming that the density distribution of the sediments in vertical direction is known based on some other a priori information such as well logging. To solve this specific problem, the geometric inversion method to recover the depth to basement has been proposed in recent years. These types of methods are mostly based on the discretization of the sedimentary columns or the sediment-basement interface instead of the 3D subsurface. By adopting this new discretization, the parameters of unknown in the inversion and the model uncertainty can be reduced significantly. In the last section of this chapter, we will discuss the most updated geometric modeling and inversion method.
2. Three-dimensional gravity modeling 2.1. Basic theory for gravity method In this section, we will introduce the basic gravity theory and the corresponding solution for solving the gravity modeling problem. We start from the basic Poisson's equation for gravity problem and introduce the solution for gravity field data caused by some mass. For numerical solution, we show the discretized formulae for calculating the gravity field data. We also briefly discuss the fast modeling methods which is based on fast Fourier transform (FFT). In certain application, the gravity problem can be reduced to study the response caused by a sediment-basement interface model (or it can be called as the density contrast model). In this case, the forward modeling problem can be solved using the column discretization or the surface discretization method based on 3D analog of Cauchy-type integral.
We first consider a simple gravity problem as shown in Figure 1 which contains a domain filled by some density distribution of ρ(r). It can be shown that the gravity field caused by this density distribution satisfies the following equation [1][2][3]: where the parameter γ denotes the gravitational constant. We can see from Eq. (1) that the gravity field g is curl free which is also referred as a conservative vector field. Take into consideration of the fact that ∇ Â (∇U) = 0, one can construct the solution for Eq. (1) by introducing the scalar gravity potential U such that: As a result, the gravity field satisfies the Poisson's equation [1,2]:

Integral representation of gravity anomaly
For geophysical problem, the observation point is usually outside the mass, in this case, Eq. (3) can be reduced to the well-known Laplace equation. In this case, the solution for gravity field outside the source domain can be written as: where D indicates the domain with anomalous mass. The corresponding Green's function for the gravity field caused by a point source can be written as follows: By comparing Eqs. (4) and (5), one can see that the gravity field caused by some mass with arbitrary shape can be calculated by integral the solution for the point source.
By taking the spatial derivative of Eq. (4), one can find the integral representation of the gravity gradiometry data which is defined as: g αβ ¼ g xx g xy g xz g yx g yy g yz g zx g zy g zz It can be proved that such gravity tensor is symmetric and the summation of the diagonal components equals to zero (this is because scalar gravity potential is Laplacian: ∇ 2 U = g xx + g yy + g zz = 0, outside the source region). As a result, among these nine gravity tensor components, only five of them are independent.
The integral formula in Eq. (4) and its derived formulae are the common approach for calculating the gravity anomaly caused by some excess mass. For numerical modeling and inversion of gravity data, it is convenient to write a discretized form of Eq. (4) by dividing the subsurface mass anomaly into a grid of prism cells: where ρ k denotes the density value of the kth cell; Δx, Δy, Δz denotes the dimension of the cell in x, y, z direction; g k is the kernel function which is defined as follows: Numerically, Eq. (8) can be evaluated using the Gaussian quadrature or the simplest central point integral formula which could be less accurate if the receiver is closer to the cell.

Fast Fourier transform (FFT) method for calculating gravity anomaly
In practical application, the modeling domain can be very large which may result in millions of cells. In this scenario, the calculation of gravity field directly using Eqs. (4) and (7) can be time and memory consuming.
Geophysicists have attempted, for decades, to apply the fast Fourier transform algorithm to potential field data modeling and inversion [4]. In the pioneering work of Parker and Oldenburg [5,6], it has been shown that the FFT method can be used to calculate the gravity field caused by a density model contains different non-flat layers. Nagendra et al. [7] have released a FORTRAN code for gravity modeling which is based on the method of Parker and Oldenburg [5,6].
The more advanced method for calculating gravity anomaly caused by complex geometry can be found in [8]. From the numerical perspective, the whole density anomaly domain is divided into a series of horizontal layers extends infinitely in x and y direction. For each layer, the corresponding gravity anomaly is calculated using FFT method and the results for each layer will be summed together to get the total gravity response. This approach usually requires a uniform discretization in x and y direction. Interested readers are recommended to get more detailed work in [8].

Differential equation for gravity modeling
Another alternative approach to overcome the large memory and computation requirement, which occurs in the analytical integral representation, is the differential equation method which includes finite difference, finite element and finite volume methods. These types of methods solve the Poisson's equation in Eq.
(3) directly with the homogeneous Dirichlet boundary condition. All these method finally result in the sparse system of equation for the scalar gravity potential. Farquharson and Mosher [9] implemented a 3D finite difference algorithm for solving the Poisson's equation of scalar gravity potential. The numerical study has shown its high accuracy (relative error is around 1%) comparing to the analytical method with direct integral over each cell. This method runs slower than the conventional integral method but it uses much less memory.
The finite element or finite volume method with unstructured mesh (e.g., tetrahedral element for three-dimensional case), can be used for solving this problem with much less memory consumption [10]. In this chapter, we will take the finite element method with unstructured tetrahedral mesh for example. We consider a tetrahedral element shown in Figure 2 with node indexing [11]. We assume that in each tetrahedral element, the scalar gravity potential is defined in the node. The gravity potential inside the tetrahedral element can be written as follows: where N j e x; y; z ð Þis the linear basis function for the tetrahedral element [11].
After applying finite element analysis to the Poisson's equation for scalar gravity potential, one can find a sparse system of equations as follows: where K is the finite element stiffness matrix [11], b is the source term which is related with the density distribution. After applying the proper boundary condition (e.g., the homogeneous Dirichlet boundary condition), the system of equation in Eq. (10) can be solved either with iterative method or the modern direct solvers. It has been demonstrated that that finite element method for gravity modeling can simulate complex geological structures. Furthermore, the method is advantageous to the conventional integral method in terms of both memory and computation speed for realistic large scale geological models.  massive sulfide deposits located in Labrador, Canada [10]. One can see that the mesh is locally refined nearby the area with complex geometry. It is difficult to simulate such complex model using the regular prism grid.

Advanced method for the modeling of density interface model
In oil and gas exploration, it is more interesting recover the favoring geological structure for oil/gas accumulation, using gravity method. Usually, the density of the basement rocks is relatively higher than the sediment density due to the compaction effect during the sedimentation process. In this scenario, the variation of the depth to basement can cause gravity anomaly on the earth's surface and can be recorded. Consider that the density of sediments and basement rocks are well known or well constrained by other geological information, we can use the observed gravity anomaly to estimate the depth to the crystalline basement. Such research topic has been investigated for decades [12][13][14][15][16][17][18]. Figure 4 shows an illustration of a synthetic sediment-basement interface model where the sedimentary rock and the basement rock are characterized by different density. The classic and most straightforward approach for solving such forward modeling problem is based on discretizing the sedimentary pack into a grid of vertical columns. The gravity field anomaly caused by each column can be calculated using the integral formula in Eq. (4). Such expressions can be further reduced in the special case of constant sediment density which results in a constant density contrast along the sediment-basement interface. In a special case of density contrast change exponentially or quadratically with depth, there exists other special formulation for the forward modeling. Interested readers can refer to [14] for more details.
The elegant theory of integral transform provides another powerful approach for solving the gravity forward modeling problem using 3D analogy of Cauchy-type integral [1]. We consider a density contrast model shown in Figure 5. The reference model with two layers are separated by a horizontal plane P (at z = À H 0 ), which is a density contrast interface. For the idealized two-layered model, the density above and below the plane P are two different constants. In real case, the actual density contrast interface Γ is an arbitrary surface. As a result, the domain D R , which is bounded by plane P and surface Γ can cause some gravity anomaly.
Using the integral transform approach [1], such gravity anomaly can be represented as follows: where h(x, y) is the relative elevation between the surface Γ and plane P at each horizontal location, Δρ is the density contrast between these two layer which is also the anomalous density inside domain D R , C Γ R (r 0 , h(x, y)d z ) is the 3D Cauchy-type integral on the surface of Γ R for the vector function h(x, y)d z , at the point of r 0 . Mathematically, the 3D analog of Cauchy-type integral for the vector function ϕ(r) can be defined as follows [1,19,20]: where S is some closed surface which bounds a 3D domain D, ϕ(r) is some vector function defined on the surface S. The 3D Cauchy-type integral itself is a vector function which satisfies the following equation outside the domain D which is bounded by surface S [1,19,20]: These important properties make it possible for using the 3D analog of Cauchy-type integral to construct solutions of gravity problem. Using the scalar representation of 3D Cauchy-type integral, one can write the solution of gravity field caused by the model shown in Figure 5 as follows [1,19,20]: where each of α, γ, η can be equal to x, y, z; P R is the projection of Γ R on the horizontal plane P. The four-index Δ symbol is defined as follows [1,19,20]: the parameter b γ is defined as follows: Note that the model described in Figure 5 is exactly a typical sediment-basement interface model when the horizontal plane P is on the earth's surface. By taking the spatial derivative of Eq. (14), one can find the expression of gravity gradiometry data using 3D analog of Cauchytype integral transform. Clearly, we can see that the 3D gravity modeling problem can be reduced to the surface integral over the density contrast surface using the Cauchy-type integral approach. As a result, the computational cost can be reduced significantly comparing to the direct integral method. As a matter of fact, this method works for a complex distribution of density contrast over vertical direction. In this case, it only requires the density contrast function in vertical direction to be integral [1,19,20].

Three-dimensional gravity inversion
As we know, the inversion of gravity data is a serious non-unique problem. In a general case, the model parameter is much larger than the observed data points. There exists infinity number of models which can fit the gravity data in a least-square sense [2]. Furthermore, the potential filed data does not have any depth resolution. In order to obtain the most reasonable solution with correct depth resolution, it is crucial to apply regularization and some a priori information during the inversion process. In this section, we will briefly discuss some modern 3D inversion techniques for gravity field data. First of all, we will formulate the inverse problem for the 3D volumetric inversion and introduce some method for minimizing the objective/parametric functional. Follow this, we will discuss the geometric inversion approach for the depth to basement inversion or density contrast interface inversion. Finally, we will introduce the joint inversion approach.

Conventional volumetric inversion
One of the most important applications of gravity inversion is to recover the 3D subsurface density distribution. The conventional approach for solving this problem is based on discretizing the subsurface into a grid of 3D prism cells. During the inversion, the density values for each cell will be adjusted in order to fit the observed gravity data. In order to fit the data with a reasonable density model, one can construct a parametric functional defined as follows [2]: where d is the vector for observed data, m is the vector for model parameters, m apr represents the a priori model, A is the forward modeling operator which is linear in this case, W d is the data weighting matrix, W m is the model weighting matrix, and α is the regularization parameter.
In practical application, a smaller regularization parameter can cause over fitting for the data and result in some unreasonable model. A larger value of regularization parameter will place strong penalty and enforce the model parameter be closer to the a priori model with the price of a bad data fitting. Several methods have been proposed for choosing the optimized regularization parameter. These methods include, but not limited to, the adaptive selection method and the L-curve method [2,21].
The model weighting matrix can be calculated based on integral sensitivity method for reasonable depth resolution [2]. In this case, the second term in the right hand side of Eq. (17) is equivalent to the minimum norm regularization which usually produce a smooth density distribution for 3D gravity inversion.
The minimization of the parametric functional in Eq. (17) can be solved using gradient type method such as steepest ascend method, conjugate gradient method and Newton method. In the gradient type inversion, one has to obtain the sensitivity matrix. For volumetric gravity inversion problem, the sensitivity matrix is exactly the same as the forward modeling matrix. For large scale inversion, the moving sensitivity domain approach is usually used to reduce the computation cost and memory consumption [22]. Within the framework of this approach, only the model cells within some distance from the receivers will be considered. If the cells are far away, the sensitivity is assumed to be zero. As a result, the full sensitivity matrix is reduced to a sparse matrix.
Usually, the inversion of gravity data in the original linear model space can result in unreasonable value of density distribution due to the absence of constraints. In reality, we usually know the estimated density values or their upper and lower boundaries. In this case, it is desirable to transfer the model parameter m in the original linear model space into the model parameterm in logarithmic space as follows [2,23]: However, the gradient type method for the deterministic inversion is characterized by the local minimum problem, which means that the solution is dependent on the selection of initial model [24]. If the initial model is not chosen properly or closer to the true model, the inverted model could be some local minimum instead of the global minimum that we usually expected.
In order to solve this problem, the methods such as Monte Carlo method, genetic algorithm, simulated annealing method can be used to reach the global minimum [24]. Comparing to the gradient type method, these methods search the optimized model from the whole model space for the global minimum. However, these methods are usually time consuming, especially for full 3D inversion which may contain millions of model parameters.
Instead of the deterministic inversion, the stochastic inversion approach can be applied to the potential field data [24][25][26]. Comparing to the deterministic inversion, the stochastic inversion can also provide uncertainty estimation for the model confidence. Recently, we start to see more publications in geophysical stochastic inversion for potential field data.

Density contrast interface inversion
The conventional prism based inversion usually produce diffusive images even though some techniques such as focusing inversion [27] can be used to enforce a sharp boundary between different geological units. In petroleum reconnaissance using gravity method, it is of great importance to estimate the depth to basement. However, from the interpretation of conventional prism based inversion, it is difficult to pick up the correct location of such interface due to the non-uniqueness of the inversion and the low resolution of inverted density distribution. In such environment, the density contrast between sedimentary rocks and basement is usually well known based on other information such as drilling. The gravity anomaly can be attributed to the variation of the sediment-basement interface. Based on the method introduced in Section 2, this type of models can be simulated using the column discretization or the Cauchy-type integral approach. In this subsection, we will mainly focus on the inversion of sedimentbasement interface using the 3D Cauchy-type integral approach.
Within the framework of this approach, we formulate the inversion with respect to the depth to basement and the density contrast value (may not necessary be a constant). Similar to the prism based inversion, we can formulate the inverse problem using the parametric functional introduced in Eq. (17). However, the model parameters now become as follows [19]: where h is the depth to basement at different horizontal locations, and ρ is the density contrast value between the sediment and basement. Comparing to the prism based inversion for the density distribution, the forward modeling operator A which is related with the 3D analog of Cauchy-type integral, for density contrast surface inversion is a nonlinear operator since the gravity data does not have a simple relationship with the depth. Such inversion can be solved with the gradient type inversion methods. The sensitivity matrix for the depth to basement and the density contrast values can be calculated by directly differentiate Eq. (14) with respect to the depth to basement and the density contrast value. In some applications, the density contrast values are usually well known based on well logging. Under this circumstance, only the depth to basement value need to be inverted and the non-uniqueness of the inversion can be reduced significantly.
As an example, we consider a 3D sedimentary interface model with the vertical section shown in Figure 6. We consider a constant density value for the basement rock. To be realistic, we assume that the density value for sediment increase exponentially with depth due to compaction. As a result, the density contrast value will decrease exponentially with depth and approach the basement density which is assumed to be a constant. The synthetic gravity data is simulated on the earth's surface and will be used for inversion. For this model, we consider that the density contrast profile with depth is already known and only invert for the depth to basement. Figure 7 shows a comparison between the true model and the inverted model using Figure 6. A vertical section (y = 0) of the sediment-interface model with exponential density contrast in vertical direction [19]. the Cauchy-type integral approach. One can see that the depth to basement and the shape of the sedimentary basin is well recovered in the inversion.
During the inversion, we use a flat surface as the initial model and a priori model. Actually, such inversion method is robust enough and does not depend too much on the selection of initial model and a priori model. However, one can use some other initial model in order to speed up the convergence of the inversion. The famous Bouguer slab formula [28] can be used as an initial model: where g B denotes the observed Bouguer gravity anomaly, Δρ 0 is the density contrast value on the earth's surface and a is the gradient of density contrast in vertical direction. As one may note, this equation works properly for the constant density contrast or linear density contrast but it does provide a good approximation of the depth to basement for a general case, e.g., exponentially increased density contrast profile with depth.

New developments in gravity inversion
In this subsection, we will briefly introduce some other new developed techniques for gravity inversion. These new methods include, but not limited to, the binary inversion, multinary inversion and the joint inversion approach.

Binary and multinary inversion of gravity data
In some application of gravity imaging such as subsurface tunnel detection, salt structure imaging, the density range of the target and the density of host rock are usually well known or well constrained [29]. However, the conventional inversion in this case will still produce a diffusive image with spread density ranges and continuous model space. In reality, the inverted density value should be clustered nearby the density value of the host rock and the density of the target.
In the case of a model with two distinct density values, the continuous inversion parameters m in the original model space can be transformed into a new binary space for inversion [30]. Zhdanov and Cox [29] introduced a multinary inversion approach for geological models with more than two density values for different geological units. Several different functions, such as Heaviside function and Gaussian function, can be used for multinary transformation [31]. Within the framework of multinary approach, the value of the inverted model parameters is enforced to be nearby the vicinities of the preselected values.

Joint inversion approach
As we know that each geophysical method and data is only sensitive to specific model parameter. Due to the inherent non-uniqueness of geophysical data inversion, the recovered model parameter from individual inversion can be ambiguous. Such ambiguity can be reduced by incorporating more a priori information into the inversion.
Different geophysical models can be coupled with each other either directly or structurally. For example, the density and seismic velocity can be related with each other by empirical equations. Alternatively, a density model can be related with the velocity model by assuming the structurally similarity. As a result, it is possible to enforce the coupling between different model parameters by inverting this different geophysical data set simultaneously. The structurally similarity based joint inversion can be achieved by minimizing the cross gradients between different model parameters [32][33][34]. Within the framework of this approach, the structural similarity, between model parameter m 1 and m 2 , can be measured by the cross gradient which is defined as follows [32]: Such regularization term will be minimized during the minimization of data misfit functional for the joint inversion approach.
Zhdanov [31] proposed a new and more flexible joint inversion approach based on Gramian constraint. Within the framework of this approach, different model parameters are coupled through the Gramian matrix which can either force the direct relationship between different model parameters or their spatial gradients. One good property of such joint inversion is that the algorithm will only enforce such coupling when it does exist and will not introduce artificial coupling when there is no relationship or coupling between different model parameters [31].
It is straightforward that the joint inversion formulation can be simplified if there exist a shared model parameter for different geophysical data set. For example, the DC electric method and the magnetotelluric (MT) method both invert for the electric conductivity. As a result, it is unnecessary to formulate the joint inversion using the approach that we have just discussed, in this scenario. Based on this idea, the joint inversion for depth to basement using different geophysical data can be greatly simplified by considering that the depth to basement is a shared model parameter for different data set such as gravity and MT data. In the meantime, each method may also have a private model parameter such as density contrast for gravity data and conductivity contrast for MT data.
Here we consider a synthetic sediment-basement interface model [19,35]. The synthetic gravity and MT data are simulated on the earth's surface for this model. These data will be used to recover the sediment-basement interface. Furthermore, we assume that the density contrast value and conductivity contrast value are also unknown. Under this circumstance, the inversion of individual data set is characterized by strong non-uniqueness. By assuming the shared model parameter of depth to basement, for gravity and MT data, during the joint inversion approach, the recovered model parameters is much closer to their true value comparing to the individual data inversion (as one can see in Figure 8).

Summary
In this chapter, we have reviewed the 3D gravity forward modeling and inversion problem. We start from the Poisson's equation for scalar gravity potential and introduce the formula of gravity field in integral form based on the Green's function which corresponds to the solution of a point source. We also introduced some other techniques, such as FFT method and differential equation method, for fast and efficient solving of the gravity forward modeling problem.
In the application of sedimentary basin modeling, we have introduced the column discretization method and the advanced Cauchy-type integral method. In the second section of this chapter, we introduced the gravity inversion problem and first start from the conventional prism based inversion. Following this, we introduce another application of gravity inversion for locating the density contrast interface. We mainly focus on the 3D inversion Figure 8. A comparison between separate gravity and MT inversion for depth to basement, and the joint inversion result for a synthetic sediment-basement interface model [19].
based on Cauchy-type integral for recovering the depth information. Finally, we have introduced some other new recently developed techniques for gravity inversion which includes, but not limits to, the multinary inversion and joint inversion approach. The interested readers are recommended to read the relevant publications cited in this chapter, for more details.