A new numerical tool for the 3D forward modeling of potential field geophysical data in the presence of rugged topography using a numerical integration scheme

Potential field geophysical data are frequently used to image geological features in volcanic systems/areas (faults, lithological contacts, alteration zones, geothermal systems, magmatic reservoirs). However, although crucial, it can prove challenging to accurately simulate data in such regions due to the major influence of strong topographic variations. To accurately account for topography with reasonable computational cost, we develop a numerical tool for the modeling and inversion of these data. The method consists of a numerical integration scheme of the integral equations predicting gravity and magnetic data on deformable hexahedral elements. The integrals are evaluated using high-order Gaussian quadrature. Physical properties of the subsurface are defined on discrete grid points, allowing to model discontinuities in the parameters not only at the surface, but also along surfaces within the models, enabling to represent faults, lithological contacts or cavities. Our method uses non-conformal meshes with automatic local refinements in regions with rapidly varying surface topography and in the vicinity of measurement points. In particular, we have developed a local and self-adaptive iterative refinement scheme based on a local convergence criterion of the numerical integration, allowing to reduce the effect of solution singularities close to observation points. The accuracy of our method is tested by comparing our model predictions with results obtained from the tomofast code (https://doi.org/10.5194/gmd-17-2325-2024) using a fine reference discretization of the topography with rectangular prisms. These tests were performed for the modeling of the gravity and magnetic effects of topography over the geothermal system of Krafla, Iceland for the case of ground-based and airborne data. Our modeling tool will ultimately be used for the independent or joint inversion of potential field data to make use of their different sensitivities in terms of physical parameters and also lateral and depth resolutions.