Three-dimensional modelling of gravity data using finite differences

Three-dimensional minimum-structure gravity inversion programs typically use a fine rectangular mesh to discretize the Earth's subsurface. The obvious forward solver for this type of Earth model is one that calculates the gravity contribution from each cell in the mesh using one of the closed-form formulae for a right rectangular prism. However, such a solver is computationally expensive, especially in memory requirements, when implemented without modification in a Newton, or Gauss–Newton, inversion procedure. Here, an alternative method is presented that involves a finite-difference solution of Poisson's equation. First, the finite-difference equations are derived for a rectangular mesh using a finite-volume approach. A linear system of equations is generated which is sparse. The solution to the system, which is the approximation to the gravitational potential, is obtained using a conjugate-gradient procedure. This maintains the sparseness of the system, and is thus efficient in memory usage. Finally, components of the gravitational acceleration are computed from the approximate potential using a finite-difference approximation of the gradient operator. The capabilities of this forward solver are illustrated by a simple example. The method is particularly suited to gradient-based inversion programs that compute the model update using an iterative procedure and evaluate the requisite matrix–vector products by means of a pseudo-forward solution.