Fast finite‐element calculation of gravity anomaly in complex geological regions

SUMMARY Forward computation of the gravity anomaly of a density model is often a necessary step in modelling the subsurface density of a region. For geologically complex regions, this step may be computationally demanding and become the bottleneck in gravity inversion. We present a fast finite-element method (FFEM) for solving boundary value problems of the gravitational field in forward computation of gravity anomaly in complex geological regions. Testing against analytical solutions show that the method is more accurate than the classical integration method in cases where density in the material body is highly heterogeneous. At the same time, FFEM is more efficient than the integration method by a factor between s/100 and s/10, where s is the number of stations at which gravity anomalies are computed. Since s is usually much greater than 10‐100 in 3-D gravity inversion in geologically complex regions, FFEM may be significantly more efficient in gravity modelling of such regions. We illustrate the utility of this method by calculating the gravity anomalies in central Taiwan. Gravity anomaly of a region is an important geophysical data set for the investigation of the subsurface density. In the search of appropriate density models, various techniques of geophysical inversion may be employed; yet forward computation of the gravity anomalies of density models is always a necessary component in the inversion process. With the advent of 3-D tomographic seismology, great complexity and heterogeneity of the subsurface geology have been revealed. The forward computation of gravity anomalies at a large number of gravity stations, as usually required in 3-D gravity inversion for the subsurface density in such complex and heterogeneous regions, may be computationally demanding even for today’s high-speed computers, and may thus become the bottleneck in gravity inversion. The classical approach to calculating the gravitational attraction of a given density distribution is by integration based on Newton’s law. In general, the gravitational attraction at an arbitrary point (x, y, z) due to density distribution ρ(ξ , η, ζ )i nav olume V may be calculated