Fast and accurate numerical modeling of gravity and magnetic anomalies is the basis of field-data inversion and quantitative interpretation. In gravity and magnetic prospecting, the computation and memory requirements of practical modeling is still a significant issue, which leads to the difficulty of using efficient and detailed inversions for large-scale complex models. A new 3D numerical modeling method for gravity and magnetic anomaly in a mixed space-wavenumber domain is proposed to mitigate the difficulties. By performing a 2D Fourier transform along two horizontal directions, 3D partial differential equations governing gravity and magnetic potentials in the spatial domain are transformed into a group of independent 1D differential equations wrapped with different wavenumbers. Importantly, the computation and memory requirements of modeling are greatly reduced by this method. A modeling example with 4,040,100 observations can be finished in approximately 28 s on a desktop using a single core, and the independent differential equations are highly parallel among different wavenumbers. The method preserves the vertical component in the space domain, and thus a mesh for modeling can be finer at a shallower depth and coarser at a deeper depth. In general, the new method takes into account the calculation accuracy and the efficiency. The finite-element algorithm combined with a chasing method is used to solve the transformed differential equations with different wavenumbers. In a synthetic test, a model with prism-shaped anomalies is used to verify the accuracy and efficiency of the proposed algorithm by comparing the analytical solution, our numerical solution, and a well-known numerical solution. Furthermore, we have studied the balance between computational accuracy and efficiency using a standard fast Fourier transform (FFT) method with grid expansion and the Gauss-FFT method. A model with topography is also used to explore the ability of modeling topography with our method. The results indicate that the proposed method using the Gauss-FFT method has characteristics of fast calculation speed and high accuracy.

Self-demagnetization due to strongly magnetic bodies can seriously affect the interpretation of magnetic anomalies. Conventional numerical methods often neglect the self-demagnetization effects and limit their use to low susceptibilities (χ<1 SI). We have developed a novel iterative method based on the integral equation and the Gauss-fast Fourier transform (FFT) technique for calculating the magnetic field from finite bodies of high magnetic susceptibility and arbitrary shapes. The method uses a segmented model consisting of prismatic voxels to approximate a complex target region. In each voxel, the magnetization is assumed to be constant, so that the integral equation in the spatial domain can reduce to a simple form with lots of merit in numerical calculation after the 2D Fourier transform. Moreover, a contraction operator is derived to ensure the convergence of the iterative calculation, and the Gauss-FFT technique is applied to reduce numerical errors due to edge effects. Our modeling results indicate that this new iterative scheme performs well in a wide range of magnetic susceptibilities (1–1000 SI). For lower susceptibilities (χ≤10 SI), the iterative algorithm converges rapidly and indicates very good correlation with the analytical solutions. At higher susceptibilities (10<χ≤100 SI), our method still performs well, but the numerical accuracy improves with a relatively slow speed. In the extreme case (χ=1000 SI), an acceptable result is also obtained after sufficient iterative computation. A further improvement in the numerical precision can be achieved by increasing the number of prismatic voxels, but at the same time, the computational time increases linearly with the size of the voxels.

Based on a brief review of forward algorithms for the computation of topographic gravitational and magnetic effects, including spatial, spectral and hybrid-domain algorithms working in either Cartesian or spherical coordinate systems, we introduce a new algorithm, namely the CP-FFT algorithm, for fast computation of terrain-induced gravitational and magnetic effects on arbitrary undulating surfaces. The CP-FFT algorithm, working in the hybrid spatial-spectral domain, is based on a combination of CANDECOMP/PARAFAC (CP) tensor decomposition of gravitational integral kernels and 2D Fast Fourier Transform (FFT) evaluation of discrete convolutions. By replacing the binomial expansion in classical FFT-based terrain correction algorithms using CP decomposition, convergence of the outer-zone computation can be achieved with significantly reduced inner-zone radius. Additionally, a Gaussian quadrature mass line model is introduced to accelerate the computation of the inner zone effect. We validate our algorithm by computing the gravitational potential, the gravitational vector, the gravity gradient tensor, and magnetic fields caused by densely-sampled topographic and bathymetric digital elevation models of selected mountainous areas around the globe. Both constant and variable density/magnetization models, with computation surfaces on, above and below the topography are considered. Comparisons between our new method and space-domain rigorous solutions show that with modeling errors well below existing instrumentation error levels, the calculation speed is accelerated thousands of times in all numerical tests. We release a set of open-source code written in MATLAB language to meet the needs of geodesists and geophysicists in related fields to carry out more efficiently topographic modeling in Cartesian coordinates under planar approximation.