3D finite-difference numerical simulation of the gravitational field using a preconditioned GMRES iterative solver

Abstract

With the evolution of geophysical surveys from traditional two-dimensional (2D) to three-dimensional (3D) models, the resulting large data volumes pose significant challenges to inversion, particularly when resolving large-scale 3D structures. A direct solver for solving an ill-conditioned linear system resulting from the finite-difference approximation of a boundary value problem requires more memory and time than iterative solvers. To overcome this limitation, an efficient iterative solver for 3D finite-difference approach is introduced to calculate the 3D gravitational potential and the associated gravitational field. Firstly, the boundary value problem associated with 3D gravitational potential is discretized using central finite-difference technique based on right rectangular prismatic grids. The resulting large unsymmetric sparse systems are then solved using the generalized minimal residual algorithm (GMRES) iterative solver in combination with incomplete LU factorization. Secondly, to obtain high-accuracy partial derivatives of gravitational potential, a high-degree Lagrange interpolation scheme is employed. Finally, three density models are applied to test the accuracy, reliability, and flexibility of our 3D finite-difference algorithm. All computational results demonstrate that our method provides an accurate approximation of the gravitational field and is applicable to 3D forward modeling.

Keywords

3D gravitational potential / boundary value problem / numerical simulation / finite-difference algorithm / GMRES iterative solver

References

[1] Hiramatsu Y, Sawada A, Kobayashi W, et al. . Gravity gradient tensor analysis to an active fault: A case study at the Togi-gawa Nangan fault, Noto Peninsula, central Japan [J]. Earth, Planets and Space. 2019, 71(1): 107

[2] Beiki M, Pedersen L B. Eigenvector analysis of gravity gradient tensor to locate geologic bodies [J]. Geophysics. 2010, 75(6): I37-I49

[3] Li W-y, Liu Y-x, Li B, et al. . Hydrocarbon exploration in the South Yellow Sea based on airborne gravity, China [J]. Journal of Earth Science. 2016, 27(4): 686-698

[4] Teknik V, Ghods A, Inan S, et al. . Potential field imaging of salt and basement structures in the southern Zagros foreland basin [J]. Pure and Applied Geophysics. 2024, 181(10): 3071-3101

[5] Mosher C R W, Farquharson C G. Minimum-structure borehole gravity inversion for mineral exploration: A synthetic modeling study [J]. Geophysics. 2013, 78(2): G25-G39

[6] Yang Y-s, Li Y-Y. Ore-controlling structures of the Qingchengzi Pb-Zn-Au-Ag orefield, northeastern China and significance for deep ore prospecting: Revealed from gravity and magnetic anomalies [J]. Ore Geology Reviews. 2023, 156: 105376

[7] Xu M L, Yang Y B, Deng Y M, et al. . Geophysics indicator of sandstone-type uranium mineralization in the northern Ordos Basin, China: Analysis from gravity and magnetic data [J]. Earth and Space Science. 2025, 12(3): e2024EA003521

[8] Fedi M, Ferranti L, Florio G, et al. . Understanding the structural setting in the Southern Apennines (Italy): Insight from gravity gradient tensor [J]. Tectonophysics. 2005, 397(12): 21-36

[9] Zhao J-j, Chen Q, Yang Y-h, et al. . InSAR, gravity, and geomagnetic observations suggest the 2020 Mw 6.5 Stanley earthquake occurred on a blind strike-slip and a normal fault [J]. Bulletin of the Seismological Society of America. 2024, 114(4): 1857-1869

[10] Pedersen L B, Bastani M. Estimating rock-vector magnetization from coincident measurements of magnetic field and gravity gradient tensor [J]. Geophysics. 2016, 81(3): B55-B64

[11] Rim H, Li Y-guo. Single-hole imaging using borehole gravity gradiometry [J]. Geophysics. 2012, 77(5): G67-G76

[12] Liu S, Hu X-y, Fedi M, et al. . Petrophysical and geophysical constrained inversion of gravity data based on starting and referenced models [J]. Journal of Geophysics and Engineering. 2025, 22(1): 36-47

[13] Roy A, Wu L-Y. Generalized Gauss-FFT 3D forward gravity modeling for irregular topographic mass having any 3D variable density contrast [J]. Computers & Geosciences. 2023, 172: 105297

[14] Li Y-g, Oldenburg D W. 3-D inversion of gravity data [J]. Geophysics. 1998, 63(1): 109-119

[15] Won I J, Bevis M. Computing the gravitational and magnetic anomalies due to a polygon; algorithms and Fortran subroutines [J]. Geophysics. 1987, 52(2): 232-238

[16] Banerjee B, Das Gupta S P. Gravitational attraction of a rectangular parallelepiped [J]. Geophysics. 1977, 42(5): 1053-1055

[17] Johnson L R, Litehiser J J. A method for computing the gravitational attraction of three-dimensional bodies in a spherical or ellipsoidal Earth [J]. Journal of Geophysical Research (1896–1977). 1972, 77(35): 6999-7009

[18] Chakravarthi V, Raghuram H M, Singh S B. 3-D forward gravity modeling of basement interfaces above which the density contrast varies continuously with depth [J]. Computers & Geosciences. 2002, 28(1): 53-57

[19] Ren Z-y, Chen C-j, Zhong Y-y, et al. . Recursive analytical formulae of gravitational fields and gradient tensors for polyhedral bodies with polynomial density contrasts of arbitrary non-negative integer orders [J]. Surveys in Geophysics. 2020, 41(4): 695-722

[20] Wu L-Y. Efficient modelling of gravity effects due to topographic masses using the Gauss-FFT method [J]. Geophysical Journal International. 2016, 205(1): 160-178

[21] Zhao G-d, Chen B, Chen L-w, et al. . High-accuracy 3D Fourier forward modeling of gravity field based on the Gauss-FFT technique [J]. Journal of Applied Geophysics. 2018, 150: 294-303

[22] Wu L-y, Tian G. High-precision Fourier forward modeling of potential fields [J]. Geophysics. 2014, 79(5): G59-G68

[23] Pan K-j, Zhang Z-h, Hu S-g, et al. . Three-dimensional forward modelling of gravity field vector and its gradient tensor using the compact difference schemes [J]. Geophysical Journal International. 2021, 224(2): 1272-1286

[24] Zhang J, Chi-Yuen W, Shi Y-l, et al. . Three-dimensional crustal structure in central Taiwan from gravity inversion with a parallel genetic algorithm [J]. Geophysics. 2004, 69(4): 917-924

[25] Cai Y-e, Wang C Y. Fast finite-element calculation of gravity anomaly in complex geological regions [J]. Geophysical Journal International. 2005, 162(3): 696-708

[26] May D A, Knepley M G. Optimal, scalable forward models for computing gravity anomalies [J]. Geophysical Journal International. 2011, 187(1): 161-177

[27] Tong X-z, Sun Y, Huang J-w, et al. . High-accuracy gravity field and gravity gradient forward modelling based on 3D vertex-centered finite-element algorithm [J]. Journal of Central South University. 2024, 31(5): 1659-1670

[28] Jahandari H, Farquharson C G. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids [J]. Geophysics. 2013, 78(3): G69-G80

[29] Farquharson C G, Mosher C R W. Three-dimensional modelling of gravity data using finite differences [J]. Journal of Applied Geophysics. 2009, 68(3): 417-422

[30] Jahandari H, Bihlo A, Donzelli F. Forward modelling of gravity data on unstructured grids using an adaptive mimetic finite-difference method [J]. Journal of Applied Geophysics. 2021, 190: 104340

[31] Hu S-g, Pan K-j, Ren Z-y, et al. . Fast large-scale gravity modeling using a high-order compact cascadic multigrid strategy [J]. Geophysics. 2023, 88(6): G125-G140

[32] Tong C-x, Tong X-z, Wei K-Q. An efficient accelerating algorithm for 3D gravity anomalies forward modeling using discrete sine transform [J]. Journal of Applied Geophysics. 2025, 241: 105831

[33] Dai S-k, Chen Q-r, Li K, et al. . The forward modeling of 3D gravity and magnetic potential fields in space-wavenumber domains based on an integral method [J]. Geophysics. 2022, 87(3): G83-G96

[34] Gharti H N, Tromp J, Zampini S. Spectral-infinite-element simulations of gravity anomalies [J]. Geophysical Journal International. 2018, 215(2): 1098-1117

[35] Martin R, Chevrot S, Komatitsch D, et al. . A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data: Application to the western Pyrenees [J]. Geophysical Journal International. 2017, 209(1): 406-424

[36] Long J-b, Farquharson C G. Three-dimensional forward modelling of gravity data using mesh-free methods with radial basis functions and unstructured nodes [J]. Geophysical Journal International. 2019, 217(3): 1577-1601

[37] Li C-w, Xiong B, Qiang J-k, et al. . Multiple linear system techniques for 3D finite element method modeling of direct current resistivity [J]. Journal of Central South University. 2012, 19(2): 424-432

[38] Zhu Z-h, Feng Y, Yang X, et al. . An efficient parallel computing method for random vibration analysis of a three-dimensional train-track-soil coupled model using Seed-PCG algorithm [J]. Journal of Central South University. 2024, 31(1): 302-316

[39] Wu X-P. A 3-D finite-element algorithm for DC resistivity modelling using the shifted incomplete Cholesky conjugate gradient method [J]. Geophysical Journal International. 2003, 154(3): 947-956

[40] Liu J-x, Liu P-m, Tong X-Z. Three-dimensional land FD-CSEM forward modeling using edge finite-element method [J]. Journal of Central South University. 2018, 25(1): 131-140

[41] Suzuki K, Fukaya T, Iwashita T. A novel ILU preconditioning method with a block structure suitable for SIMD vectorization [J]. Journal of Computational and Applied Mathematics. 2023, 419: 114687

[42] Saad Y. Iterative methods for sparse linear systems [M]. 2003, Philadelphia, SIAM Press

[43] Baker A H, Jessup E R, Manteuffel T. A technique for accelerating the convergence of restarted GMRES [J]. SIAM Journal on Matrix Analysis and Applications. 2005, 26(4): 962-984

[44] Dubey C P, Tiwari V M. Computation of the gravity field and its gradient: Some applications [J]. Computers & Geosciences. 2016, 88: 83-96

[45] Garcia-Abdeslem J. Gravitational attraction of a rectangular prism with depth-dependent density [J]. Geophysics. 1992, 57(3): 470-473

[46] Garcia-Abdeslem J. Short Note: The gravitational attraction of a right rectangular prism with density varying with depth following a cubic polynomial [J]. Geophysics. 2005, 70(6): J39-J42

[47] Chen C-j, Ren Z-y, Pan K-j, et al. . Exact solutions of the vertical gravitational anomaly for a polyhedral prism with vertical polynomial density contrast of arbitrary orders [J]. Geophysical Journal International. 2018, 214(3): 2115-2132

[48] Chakravarthi V, Sundararajan N. 3D gravity inversion of basement relief: A depth-dependent density approach [J]. Geophysics. 2007, 72(2): I23-132

RIGHTS & PERMISSIONS

Central South University