Forward modelling of gravity and magnetic data in the presence of topography using a numerical integration scheme with deformable elements

SUMMARY Volcanic and geothermal areas usually feature highly rugged and variable topography, with both steep and relatively flat areas. Potential field geophysical data being very sensitive to topographic changes, it can be challenging to accurately simulate data in such regions. Traditional modelling approaches for potential field data typically involve discretizing the physical space using rectangular prisms, or alternatively using slightly more sophisticated geometries such as polyhedrons. However, in both cases, these discretizations may lead to an under- or overestimation of topographic effect on simulated data, unless a very fine discretization around observation points is used, consequently leading to an increase in the computational burden. To address this issue, we introduce new methodological strategies to simulate potential field data in regions with rugged topography. We propose the use of a numerical integration scheme over deformable hexahedral elements that conform to topographic variations, to numerically evaluate the integral equations governing gravity and magnetic anomalies. Density and magnetization properties are defined at grid nodes and interpolated within the elements using polynomial functions. The numerical integration relies on the use of a Gaussian quadrature, combined with a local refinement of the mesh designed to simultaneously increase the precision of the numerical integrals by reducing the effect of singularities and incorporating a more precise description of the topography at the vicinity of measurement points. This refinement is coupled with a discontinuous optimization of the elevation of geometrical nodes lying at the surface of the model by minimizing the distance between the interpolated topographic surface and a point cloud provided by a high-resolution Digital Elevation Model (DEM). The predictions at each observation points are computed with different refined versions of the mesh in order to minimize computational cost. Physical properties within the refined elements are extrapolated from the nodes of the original non-refined mesh in order not to introduce additional unknown physical parameters for the later implementation of this code for the inversion of the potential field data. The approach is first validated by comparing gravity and magnetic predictions with the analytical solution for a simple rectangular prism. We further test the method with realistic simulations using the complex topography of the Krafla geothermal area in northern Iceland. Simulation results are compared against the forward gravity and magnetic responses of a high-resolution benchmark model (discretized with 2 m $\times$ 2 m rectangular prisms) computed using the open-source software Tomofast-x. Finally, we use the same high-resolution model to quantify typical errors associated with coarser rectangular prism discretizations (cell sizes of 25, 50 and 100 m). We demonstrate that these errors can be comparable in magnitude to real geophysical anomalies and are not reproduced when using our proposed numerical approach.