Backend-agnostic Julia framework for 3D modeling and inversion of gravity data

TL;DR

3D gravity inversion is inherently ill-posed and computationally intensive due to non-uniqueness and large model spaces. The authors implement a data-space inversion framework in Julia with a backend-agnostic kernel abstraction, enabling CPU and GPU execution and reducing the problem from model space size $M$ to data space size $N$, via updates like $\delta \mathbf{m} = \mathbf{W}_m \mathbf{C}^\top (\mathbf{C} \mathbf{W}_m \mathbf{C}^\top + \mathbf{W}_d)^{-1} \delta \mathbf{d}$. Forward modeling employs a vertex-centered prism formulation for gravity, while inversion minimizes a regularized objective and leverages a matrix identity to avoid inverting the full $M\times M$ system. GPU acceleration yields substantial speedups (up to 1–2 orders of magnitude) for large models (up to ≈$3.3\times10^6$ prisms), and synthetic field tests demonstrate recovery of complex structures (vertical/dipping dykes) with geologically meaningful constraints. The framework is validated on field data, showing strong concordance with independent geological information and highlighting the potential for scalable, extensible, high-resolution gravity inversion workflows in HPC environments.

Abstract

This paper presents a high-performance framework for three-dimensional gravity modeling and inversion implemented in Julia, addressing key challenges in geophysical modeling such as computational complexity, ill-posedness, and the non-uniqueness inherent to gravity inversion. The framework adopts a data-space inversion formulation to reduce the dimensionality of the problem, leading to significantly lower memory requirements and improved computational efficiency while maintaining inversion accuracy. Forward modeling and inversion operators are implemented within a backend-agnostic kernel abstraction, enabling execution on both multicore CPUs and GPU accelerators from a single code base. Performance analyses conducted on NVIDIA CUDA GPUs demonstrate substantial reductions in runtime relative to CPU execution, particularly for large-scale datasets involving up to approximately 3.3 million rectangular prisms, highlighting the scalability of the proposed approach. The inversion incorporates implicit model constraints through the data-space formulation and depth-weighted sensitivity, which mitigate depth-related amplitude decay and yield geologically coherent, high-resolution subsurface density models. Validation using synthetic models confirms the ability of the framework to accurately reconstruct complex subsurface structures such as vertical and dipping dykes. Application to field gravity data further demonstrates the robustness and practical utility of the GPU-accelerated framework, with the recovered models showing strong consistency with independent geological constraints and prior interpretations. Overall, this work underscores the potential of GPU-enabled computing in Julia to transform large-scale gravity inversion workflows, providing an efficient, extensible, and accurate computational solution for high-resolution geophysical studies.

Figures (12)

• Figure 1: Computational workflow of the 3D gravity inversion framework
• Figure 2: Log-log plot comparing computation time for sensitivity matrix construction using GPU versus CPU implementations across varying total number of cells. The GPU demonstrates dramatically better scaling, achieving significant speedup for large-scale models up to $3.3 \times 10^{6}$ cells, making GPU acceleration essential for practical 3D gravity inversion workflows.
• Figure 3: Combined visualization of the synthetic model and inversion results. (a) Synthetic density model displayed using an isosurface representation. (b) Inverted density model visualized with isosurfaces for comparison with the true geometry. (c) Slice view of the synthetic model showing the internal density distribution. (d) Slice view of the inverted model illustrating the recovered subsurface structure. Together, the panels provide a direct comparison between the true density distribution and the inversion output, demonstrating the ability of the method to recover both the geometry and internal characteristics of the anomalous bodies.
• Figure 4: Combined visualization of the synthetic model and inversion results for the dipping dyke example. (a) Isosurface representation of the synthetic density model. (b) Isosurface visualization of the inverted model showing the recovered dyke geometry. (c) Slice view of the synthetic model illustrating internal density structure. (d) Slice view of the reconstructed density model obtained using data space inversion. Together, the panels provide a comprehensive comparison between the true model and the inversion output, demonstrating the ability of the method to resolve dipping subsurface structures.
• Figure 5: Geological map of the Tangarparha study area, Odisha, India, illustrating the distribution of lithological formations—such as sheared granite, quartzofeldspathic gneiss, and mafic/ultramafic rocks—and key structural features.