Robust 3D Joint Inversion of Gravity and Magnetic Data: A High-Performance Computing Approach

: One of the fundamental challenges in geophysics is the calculation of distribution models for physical properties in the subsurface that accurately reproduce the measurements obtained in the survey and are geologically plausible in the context of the study area. This is known as inverse modeling. Performing a 3D joint inversion of multimodal geophysical data is a computationally intensive task. Additionally, since it involves a modeling process, finding a solution that matches the desired characteristics requires iterative calculations, which can take days or even weeks to obtain final results. In this paper, we propose a robust numerical solution for 3D joint inversion of gravimetric and magnetic data with Gramian-based structural similarity and structural direction constraints using parallelization as a high-performance computing technique, which allows us to significantly reduce the total processing time based on the available Random-Access Memory (RAM) and Video Random-Access Memory (VRAM)and improve the efficiency of interpretation. The solution is implemented in the high-level programming languages Fortran and Compute Unified Device Architecture (CUDA) Fortran, capable of optimal resource management while being straightforward to implement. Through the analysis of performance and computational costs of serial, parallel, and hybrid implementations, we conclude that as the inversion domain expands, the processing speed could increase from 4× up to 100× times faster, rendering it particularly advantageous for applications in larger domains. We tested our algorithm with two synthetic data sets and field data, showing better results than standard separate inversion. The proposed method will be useful for joint geological and geophysical interpretation of gravimetric and magnetic data used in exploration geophysics for example minerals, ore, and petroleum search and prospecting. Its application will significantly increase the reliability of physical-geological models and accelerate the process of data processing.


Introduction
Geophysical methods are one of the main ways to study the Earth's crust and search and explore the mineral resources in the Earth's interior [1][2][3].Geophysical methods can be classified based on the physical fields used in their application.Gravity and magnetic methods thus belong to the so-called potential methods of geophysics since they are based on potential theory [4][5][6][7].
One of the fundamental challenges in geophysics is the calculation of distribution models of physical properties in the subsurface that can accurately reproduce measurements obtained in the field, in turn, these models are geologically possible in the context of the study area, which is known as inverse modeling.Any inversion problem is subject to three constraints: (1) the solution must exist, (2) the process is unstable, and (3) the solution is not unique [4,[8][9][10][11].
The informativeness and reliability of the interpretation of geophysical data depend greatly on the method of their inversion.Many studies have been devoted to this topic; Appl.Sci.2023, 13, 11292 2 of 24 we will briefly try to summarize different techniques of joint inversion, their advantages, and their disadvantages.At the same time, it is important to realize that the inversion of different methods may differ in the mathematical apparatus used for their joint processing.
For example, in petrophysical inversion [12,13], a reliable result can be achieved using Newton-Gauss iteration [14].The authors note that separate inversion of electromagnetic or seismic data is highly ambiguous and leads to great ambiguity in determining porosity and fluid saturation.Using the Gauss-Newton algorithm equipped with multiplicative regularization and a proper data-weighting scheme in the example with real data, the authors showed that up to 20% errors in the saturation and porosity exponents in Archie's equations do not cause significant errors in the inversion results [12,13].
Spichak [15] in his work on this topic points out that in petrophysical data inversion, it is necessary to apply two levels of initial data, the accuracy and representativeness of which influence the final result.Lelievre [16] introduced so-called weight coefficients (regularization terms) for this purpose.As a result, when applying this type of joint inversion, it is necessary to introduce two regularization levels.Level 1 (conductivity, seismic velocity, and density) and Level 2 (porosity and water content).The disadvantage of this approach is that a lack of accuracy or number of measurements of Level 1 parameters will result in an unreliable inversion as a whole.
Franz et al. [17] analyze the compatibility of different methods in the joint inversion of geophysical data in a study of the Namibian continental margin.They note that electrical resistivity models obtained by inversion of marine data at shallow water depths are characterized by poor resolution for lithologic boundaries.To address this important problem, the authors use a joint inversion of magnetotelluric, gravimetric (fixed structural density model coupling with satellite gravity data), and seismic (fixed gradient velocity model) data.The number of iterations varies from 283 to 380; the Final MT and gravity data RMS are 3.01 and 1.52, respectively.
In turn, Gallardo & Meju [18] proposed a joint inversion of magnetotelluric and seismic data, applying this method to the interpretation of geologic features such as sedimentary basins, crustal and upper mantle structures, oil reservoirs, etc.The joint inversion of these data with added noise has led to good results, but its numerical implementation is often unstable and requires regularization [15].The main disadvantage of this method is the need for structural semblance of the models and too much computer core memory, especially for three-dimensional inversion.
It is important to note that inversion can be interactive or sequential.Dell'Aversana [19,20] describes this process in detail, using the construction of a velocity model following a tomographic inversion.If necessary, the new seismic section is further used to simulate resistivity, etc.The approach is advantageous as it is relatively simple and requires much less computer core memory than the simultaneous joint inversion [15].
According to Reimann et al. [21], inversion methods can be categorized as follows: Cluster analysis is grouping according to the proximity of samples in the space of properties [22,23]; the Gaussian classifier is based on the assumption that the vector of properties has an n-dimensional Gaussian distribution of the conditional pdf within each lithological group [23]; K-means clustering implies an iterative search of a set of points (centers) so as to minimize the mean squared distance from each data point to its nearest center in training samples [24]; Discriminant analysis searches the functions of properties that separate groups in an optimal way (in terms of root mean square) in order to choose the combination of linear coefficients to maximize the variance of group centroids while minimizing the variance within the groups [25].
In this regard, the approach of the objective function according to the type of problem to be solved and the way to obtain a solution are fundamental to achieving satisfactory results; thus, a mathematical approach to obtain a solution capable of fitting the needs of each particular case is what we call "robust".The evaluation of uncertainty is an important step in inversion [26] since inverse problems always belong to decision-making processes and, by nature, are ill-posed, i.e., there are different solutions (called equivalents) that are compatible with the prior information and fit the observed data within the same error bounds [11,[27][28][29].
The main objective of geophysical data inversion, and in this case, gravity and magnetic data, is to solve the inverse problem of geophysics and related geological problems.Such problems include geologic structure, tectonics and geodynamics, volcanism and seismic processes, hydrogeology and mineral resource studies, etc.The advantages of these methods are relatively low cost, high productivity of work, and the possibility to obtain information about large territories in great depth.Many examples of successful applications of potential methods in different areas of the world and for solving both prospecting and theoretical geologic problems are described [30][31][32].In the first stages, the inversion of geophysical data is carried out in two dimensions and for each method separately [33][34][35][36][37][38].The next step may be three-dimensional modeling [39,40].The last step may be a joint inversion by two or more methods, performed in two and three dimensions [41][42][43].A serious disadvantage of gravity and geomagnetic methods is the ambiguity of the interpretation of the obtained data.Beyond the reductions as such, data analysis, functional fitting, and regional-local anomaly separation will also introduce errors; however, these remain essentially undefined because the basic assumptions are arbitrary and regional-local separation as such cannot be made in an unambiguous manner [44][45][46].It should be noted that this is a problem for most geophysical methods.One way to overcome this problem, or at least to reduce uncertainty about the results, is to apply different filters and integrate methods [47,48].However, this is not always successful either.The best results can currently be achieved by the joint inversion of several methods, significantly improving the reliability of the interpretation [49][50][51][52][53]. On the other hand, it makes the mathematical procedure of data processing more complex and leads to the need for higher time consumption and more computational power; high-performance computing techniques are an important task to improve the mathematical process of joint data processing in order to obtain solutions in a reasonable time with ordinary personal computers.
High-performance computing (HPC) techniques to obtain solutions to inversion problems have been used since the late 1980s, as mentioned before [54].However, starting in the second decade of the 2000s, a new era of HPC has developed thanks to the arrival of the new graphic processing units (GPU), its high parallelization capacity, and easy access outside the business or academic field due to its low cost.The most commonly used techniques in geophysics have focused mainly on electromagnetic and seismic methods and can be generally classified as parallelization of numerical methods [55][56][57][58][59][60][61], based on neural networks and deep learning [62][63][64][65], and more recently, based on quantum computing [66].
In this case, the variational modeling approach has been employed to construct an objective function or parametric functional consisting of (1) a mismatch term measuring the distance between F(x) and y, usually via the L2 norm, and (2) a regularization functional favoring appropriate minimizers and penalizing potential solutions with undesirable structures [67] using the form: The regularization functional R(x) represents a family of linear operators that should stabilize the problem and bound desirable solutions within the infinite space of possible solutions, while α is a scalar called the regularization parameter that gives weight to the regularization functional.The higher this value, the more influence it has on the solution; conversely, the closer to zero this value is, the closer it is to a solution without regularization.The purpose of regularization tools is to provide the user with easy-to-use routines based on numerically robust and efficient algorithms for doing experiments with regularization of discrete, ill-posed problems [68].Two widely used regularizers are the Laplacian, which adds smoothness to the models because it is known that (general) changes in properties occur gradually; therefore, smoothed distributions are expected.On the other hand, sometimes there is subsurface information obtained directly or indirectly in previous studies, and the solution should incorporate such information, so an a priori model is used as a second regularizer.In addition, it is possible to know the predominant orientations of structures in the study area and some directional preference regularizer may be useful; therefore, the implementation of such regularizers allows robust solutions capable of fitting a wide variety of needs according to geological-structural context.
Joint inversion as a regularizer to mitigate the non-uniqueness of solutions has proven to be very effective [69][70][71][72][73][74][75][76].This assumes that information obtained from different methods, whether directly or indirectly, comes from the same source interacting with different physical fields.Therefore, it is reasonable to expect inversion results where the distribution of different physical properties in the subsurface of the same study area have some relationship, as well as to propose the generation of models that feedback through each other through correspondence relationships that can be (1) intrinsic to rock physics or (2) structural.Zhdanov et al. [73] shows that in a Gramian space, the norm is a function that provides a measurement of the correlation between that function and the different parameters of the model within a Gram matrix.Due to its quadratic form similar to the L2 norm, it can be used as a regularizer of structural similarity, with the advantage that this function is differentiable, thus avoiding resorting to some linearization as in the case of crossed gradients [41].
To perform a 3D inversion of multiple types of geophysical data together is a task with an extremely high computational cost that grows exponentially as the resolution of the solutions increases.Furthermore, since it is a modeling process, finding a solution that meets the desired characteristics involves calculating solutions iteratively, which can take several days or weeks to obtain a final result.
In this work, we propose a robust numerical solution for 3D joint inversion for gravimetric and magnetic data using Gramian-based structural similarity and greater structural direction constraints, using parallelization as a high-performance computing technique, which allows us to significantly increase the processing speed and improve the efficiency of interpretation.

Methods
This methodology describes a 3D joint inversion for gravimetric and magnetic data with smoothness, a priori information, directional preference, and Gramian-based structural similarity constraints.
G ∇m (2) , ∇m (1)     ; where I is an identity matrix, A is the linear operator for the direct modeling of gravimetry and magnetometry respectively, m is the vector encoding the model parameters or physical property distribution of a 3D domain (subsurface), d 0 is the vector containing the measured data in a 2D domain (surface), m apr is the vector containing a priori information about the distribution of properties in the subsurface, W d is the weights diagonal matrix of inverse standard deviation assigned to each data and W m is the weights diagonal matrix of inverse standard deviation assigned to each cell in the model.The regularization parameters α (1) , α (2) , β (1) , β (2) , γ (1) , γ (2) , η (1) , η (2) , µ (1) , µ (2) determine the weight of the regularization functionals, which allows varying importance to be given to the regularizers in the inversion solutions for each specific case, thus having more control over the geometries of the contrasts obtained in the inversion models, thereby leaving the freedom to choose what kind of results are expected to the user's objectivity.

The linear operators L, M, and N correspond to the Laplacian
which smoothes the distribution of properties in the solutions, , which is the directional derivative of the strike θ and pitch φ angles and whose function is to specify a preferential structural direction, as well as the vertical derivative N = ∂ f ∂z , which depending of its regularization parameter sign favors vertical or horizontal properties distributions in the subsurface.These matrices were computed from the linearization of the 3D domain and discrete numerical derivatives using central differences and forward and backward differences depending on which position of each node in the 3D domain (Figure 1).This allowed us to obtain the curvature and derivatives of the domain without the need to consider extra cells or reduce the domain.
The operators A (1) and A (2) were calculated using analytical solutions for straight prisms following the development of [77,78] for gravimetry and magnetometry, respectively (See Appendix C).
The Gramian G [73] is defined as the determinant of the Gram matrix for the elements m (1) and m (2) .The equivalence of the Gramian with respect to the L2 norm of crossed gradients is detailed in Appendix B of [76].

Gravity and Magnetic Inversion Using Conjugate Gradient Algorithm
The following functional was proposed to be minimized through an iterative process to approximate the solution using the conjugate gradient method: Thus, summarizing the iterative scheme is defined as follows: where the stopping conditions will be (1) the maximum number of iterations and (2) the minimum misfit error achieved, both configured at the beginning of the process.

Coding Guidelines
For the computational implementation, memory management, the number of calculations to be performed, and the size of the arrays involved in operations were considered important points to optimize and achieve better performance.The Fortran 2008 standard [79] was used as the programming language because (1) it allowed the management of the memory size used by each array in a simple way, (2) also allowed to implement of the parallelization of matrix-vector type multiplications through the recently (2021) released free Compute Unified Device Architecture (CUDA) Fortran libraries within the Nvidia Corporation (NVIDIA) (Santa Clara, CA, USA) high-performance computing software development kit (HPC SDK v23.9), (3) as well as better performance compared to other high-level programming languages without complicating its parallel programming thanks to the similar syntax between the language and the parallelization library in contrast with other low-level programming languages.CUDA is a general parallel computing library developed by NVIDIA, which is a powerful tool compatible with popular programming languages.
Although it is convenient to process with concatenated arrays and vectors to describe the procedure, it was decided to perform the calculations of each type of geophysical data separately.This allowed to reduce the size of all the matrices to be multiplied by 50% since 2/4 parts of the concatenated matrices are filled with zeros.This implied performing a larger number of operations, but smaller in size, allowing faster processing and using a smaller amount of memory to store the temporary arrays.Even with these considerations, the memory usage and processing time grow exponentially as the inversion domain increases, depending on the number of cells to be processed and not on their dimensions.
In order to take full advantage of the computational power available during the modeling process, three different routines were programmed: A single-core version on the Central Processing Unit (CPU) with the amount of RAM available as the only limiting factor for the inversion domain size at the request of very long processing times, which has been compiled with two different compilers, gfortran and nvfortran, and showed consistently to be at least 4× times faster (Table 1).A parallel version that uses the Graphics processing unit (GPU) to perform all matrix-vector multiplications in the process in parallel, which is many times faster at the cost of relying on the amount of available VRAM memory, a much more difficult computational resource to obtain in large quantities.Finally, a CPU-GPU hybrid version aims to be a balance between speed and resource usage by performing only the calculations of structural similarity regularizer in parallel.Table 1 summarizes the performance of different processing routines for different inversion domain sizes based on RAM usage, VRAM usage, processing time, and how many times faster it is compared to the slower single-core version.Meanwhile, a visual comparison is presented in Figure 2. Performance tests were conducted on an easy to access desktop hardware (CPU: Ryzen 5

Results and Validation
To demonstrate the capabilities of a more robust parametric functional (Equation ( 2)), a comparison of results with respect to the standard separate inversion functional (Equation ( 17)) was performed: The following expressions [80] were used to evaluate the RMS of normalized residuals and the model convergence for each iteration for each data type: For gravimetric data when j = 1 and magnetic data when j = 2 and where C is the covariance matrix associated with the data, n is the number of records corresponding to each type of geophysical data and ε is a small positive arbitrary value to avoid division by zero.

Synthetic Test Model 1
To validate the correct operation of the parametric functional, a synthetic experiment proposed by [81,82] was replicated, consisting of a dike inclined at 45 degrees and buried at 50 m deep in a half-space of zero contrast, with constant density and magnetization of 1 g/cm 3 and 1 A/m, respectively.This model is made up of 20 × 20 × 10 (4000) voxels of 50 m per side, forming a volume of 1000 × 1000 × 500 cubic meters (Figure 3), from which gravity and magnetic anomalies were calculated by direct modeling, obtaining maps of 20 × 20 cells (400 data points) separated each 50 m and contaminated by Gaussian noise with a standard deviation of 2% to verify the sensitivity to noise of the inversion solutions.For the magnetic case, a field strength B 0 = 40, 000, I = 45 • and D = 45 • was considered.Figure 4 shows the resulting gravimetric and magnetic anomalies that were used as input data for the inverse modeling.

Synthetic Test Model 1
To validate the correct operation of the parametric functional, a synthetic experiment proposed by [81,82] was replicated, consisting of a dike inclined at 45 degrees and buried at 50 m deep in a half-space of zero contrast, with constant density and magnetization of 1 g/cm 3 and 1 A/m, respectively.This model is made up of 20 × 20 × 10 (4000) voxels of 50 m per side, forming a volume of 1000 × 1000 × 500 cubic meters (Figure 3), from which gravity and magnetic anomalies were calculated by direct modeling, obtaining maps of 20 × 20 cells (400 data points) separated each 50 m and contaminated by Gaussian noise with a standard deviation of 2% to verify the sensitivity to noise of the inversion solutions.For the magnetic case, a field strength  = 40,000, I = 45° and D = 45° was considered.Figure 4 shows the resulting gravimetric and magnetic anomalies that were used as input data for the inverse modeling.

Standard Separate versus Robust Joint Inversion Results of Test Data 1
This experiment was performed with data contaminated by Gaussian noise to consider its effect on the inversion results.We considered hypothetical clipping information from two wells with real density values as a priori information for the gravity case the Well_1 at 350-400 m and Well_2 at 650-700 m on the x-axis in order to evaluate how the magnetic case improves with joint inversion.40 iterations were run, and the stopping condition was convergence between models of less than 1%.The regularizers shared in standard separate (Equation ( 17)) and robust joint inversion functional (Equation ( 2)) were kept constant.For the joint inversion case, the "structural direction" regularizer was also used to indicate the strike and dip angles of the structure.Table 2 specifies the rest of the inversion parameters.
The results of the experiment at epoch 30 are plotted in Figure 5 with a general color scale ranging from 0 to 1 g/cm 3 and A/m for the gravimetric and magnetic cases, respectively.Figure 6 also plots the model convergence of each iteration in percentage, which represents how much a model changes with respect to the past iteration and the RMS error of normalized residuals, which reflects the mismatch between the input synthetic data and the data calculated from the inversion models of each iteration.18)) and data misfit (Equation ( 19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

We have shared access to an alpha version of the program under development at supplementary materials section where this synthetic model can be reproduced. All needed instructions are found within the repository description.3.2. Synthetic Test Model 2
A second test model with bodies in different positions and property contrast, inspired by the synthetic models of [61,83], was used to further test the operation of the proposed parametric function in more complex situations.It is made of 30 × 20 × 15 (9000)    18)) and data misfit (Equation ( 19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

We have shared access to an alpha version of the program under development at supplementary materials section where this synthetic model can be reproduced. All needed instructions are found within the repository description.3.2. Synthetic Test Model 2
A second test model with bodies in different positions and property contrast, inspired by the synthetic models of [61,83], was used to further test the operation of the proposed parametric function in more complex situations.It is made of 30 × 20 × 15 (9000)  18)) and data misfit (Equation ( 19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.
We have shared access to an alpha version of the program under development at supplementary materials section where this synthetic model can be reproduced.All needed instructions are found within the repository description.

Synthetic Test Model 2
A second test model with bodies in different positions and property contrast, inspired by the synthetic models of [61,83], was used to further test the operation of the proposed parametric function in more complex situations.It is made of 30 × 20 × 15 (9000) voxels of 50 m per side, forming a volume of 1500 × 1000 × 750 cubic meters (Figure 7), from which gravity and magnetic anomalies were calculated by forward modeling, obtaining maps of 30 × 20 cells (600 data points) separated by 50 m; also, these maps were contaminated with Gaussian noise with a standard deviation of 5%.For the magnetic case, a field strength B 0 = 40, 000, I = 45 • and D = 15 • was considered.Figure 8 shows the resulting gravimetric and magnetic anomalies that were used as input data for the inverse modeling.
Appl.Sci.2023, 13, x FOR PEER REVIEW 13 of 27 voxels of 50 m per side, forming a volume of 1500 × 1000 × 750 cubic meters (Figure 7), from which gravity and magnetic anomalies were calculated by forward modeling, obtaining maps of 30 × 20 cells (600 data points) separated by 50 m; also, these maps were contaminated with Gaussian noise with a standard deviation of 5%.For the magnetic case, a field strength  = 40,000, I = 45° and D = 15° was considered.Figure 8 shows the resulting gravimetric and magnetic anomalies that were used as input data for the inverse modeling.

Standard Separate versus Robust Joint Inversion Results of Test Data 2
For this inversion test, maps generated by forward modeling of test model 2 were used as input data automatically generating a domain of 32 × 22 × 15 (10,560) voxels of 50 m on each side.Moreover, 20 iterations were run, and the stopping condition was convergence between models less than 1%.Only a smoothness regularizer was used in standard separate inversion (Equation ( 17)), while for robust joint inversion (Equation ( 2)), more information was given to demonstrate how more information apart from a priori data can be crucial to obtain better results.The parameters used in the inversion are described in Table 3. voxels of 50 m per side, forming a volume of 1500 × 1000 × 750 cubic meters (Figure 7), from which gravity and magnetic anomalies were calculated by forward modeling, obtaining maps of 30 × 20 cells (600 data points) separated by 50 m; also, these maps were contaminated with Gaussian noise with a standard deviation of 5%.For the magnetic case, a field strength  = 40,000, I = 45° and D = 15° was considered.Figure 8 shows the resulting gravimetric and magnetic anomalies that were used as input data for the inverse modeling.

Standard Separate versus Robust Joint Inversion Results of Test Data 2
For this inversion test, maps generated by forward modeling of test model 2 were used as input data automatically generating a domain of 32 × 22 × 15 (10,560) voxels of 50 m on each side.Moreover, 20 iterations were run, and the stopping condition was convergence between models less than 1%.Only a smoothness regularizer was used in standard separate inversion (Equation ( 17)), while for robust joint inversion (Equation ( 2)), more information was given to demonstrate how more information apart from a priori data can be crucial to obtain better results.The parameters used in the inversion are described in Table 3.

Standard Separate versus Robust Joint Inversion Results of Test Data 2
For this inversion test, maps generated by forward modeling of test model 2 were used as input data automatically generating a domain of 32 × 22 × 15 (10,560) voxels of 50 m on each side.Moreover, 20 iterations were run, and the stopping condition was convergence between models less than 1%.Only a smoothness regularizer was used in standard separate inversion (Equation ( 17)), while for robust joint inversion (Equation ( 2)), more information was given to demonstrate how more information apart from a priori data can be crucial to obtain better results.The parameters used in the inversion are described in Table 3.
The results of the experiment at epochs 1, 10, and 20 are plotted in Figure 9 with a general color scale ranging from −1 to 1 g/cm 3 and A/m for the gravimetric and magnetic case, respectively.Figure 10 also plots the model convergence of each iteration in percentage, which represents how much a model changes with respect to the past iteration and the

Field Data Testing
Figure 10.Model convergence (Equation ( 18)) and data misfit (Equation ( 19)) achieved in inversion of test data 2 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

Field Data Testing
To test the program and the robust joint inversion algorithm, we performed the inversion of preliminary field data from Los Contreras phreatomagmatic maar (Figure 11), located at the San Luis Potosí volcanic field, Central Mexico [84,85].Gravity data were obtained at 50 m through walks using a CG-5 Autograv gravity meter with a reading resolution of 0.001 mGal.Measurement of the Earth's magnetic field was realized at the same spots using a GSM-19 high-sensitivity Overhauser effect magnetometer with 0.01 nT resolution and 0. To test the program and the robust joint inversion algorithm, we performed the inversion of preliminary field data from Los Contreras phreatomagmatic maar (Figure 11), located at the San Luis Potosí volcanic field, Central Mexico [84,85].Gravity data were obtained at 50 m through walks using a CG-5 Autograv gravity meter with a reading resolution of 0.001 mGal.Measurement of the Earth's magnetic field was realized at the same spots using a GSM-19 high-sensitivity Overhauser effect magnetometer with 0.01 nT resolution and 0.  After corrections were performed to the raw data, the residual gravity anomaly was calculated by subtraction of the regional components from the observed complete Bouguer anomaly, and simultaneously, the regional magnetic field was subtracted from the observed magnetic anomaly data.After data reduction, we statistically interpolate the data using variograms to obtain better coverage of the study area and resample a regular mesh at 100 m of spacing, obtaining 128 data points for each (gravity and magnetic) method (Figure 12).After corrections were performed to the raw data, the residual gravity anomaly was calculated by subtraction of the regional components from the observed complete Bouguer anomaly, and simultaneously, the regional magnetic field was subtracted from the observed magnetic anomaly data.After data reduction, we statistically interpolate the data using variograms to obtain better coverage of the study area and resample a regular mesh at 100 m of spacing, obtaining 128 data points for each (gravity and magnetic) method (Figure 12).We performed a joint inversion with a smoothness regularizer for 20 iterations with a stopping condition of convergence between models less than 1%.No a priori information, structural direction, variable standard deviation to data or models was introduced due to better results according to the area geological-structural context, which are beyond the scope of this article.A domain of 21 × 11 × 10 (2310) voxels was automatically generated based on the spatial limits of the input data and the voxel size of 100 m per side.Results are shown in Figure 13, slicing density and magnetization domains from top to bottom at 300 m until no more physical property contrast is visible.Additionally, the convergence and RMS graphs are also displayed in Figure 13.We performed a joint inversion with a smoothness regularizer for 20 iterations with a stopping condition of convergence between models less than 1%.No a priori information, structural direction, variable standard deviation to data or models was introduced due to better results according to the area geological-structural context, which are beyond the scope of this article.A domain of 21 × 11 × 10 (2310) voxels was automatically generated based on the spatial limits of the input data and the voxel size of 100 m per side.Results are shown in Figure 13, slicing density and magnetization domains from top to bottom at 300 m until no more physical property contrast is visible.Additionally, the convergence and RMS graphs are also displayed in Figure 13.

Discussion
Here we list the main discussion issues related to the solution to the problem posed in the paper.In brief, they can be summarized as follows: It is possible to further optimize the computational implementation, e.g., by CPU parallelization using MPI; however, this only decreases the processing time (objective already met) without enhancing memory management.On the other hand, by incorporating

Discussion
Here we list the main discussion issues related to the solution to the problem posed in the paper.In brief, they can be summarized as follows: It is possible to further optimize the computational implementation, e.g., by CPU parallelization using MPI; however, this only decreases the processing time (objective already met) without enhancing memory management.On the other hand, by incorporating variable resolution both horizontally and vertically through cells of varying dimensions, we would be able to differentiate the domain more efficiently, thereby reducing memory usage.
At this point in development, the program is only capable of working with an equal amount of data for both types of input data; likewise, the resolution of the joint inversion models is the same for both geophysical techniques, limiting its application to the spatial coverage of both data types.Additionally, the way we construct our 3D domain does not consider topography yet.
The shown processing times are representative, as they depend on the quantity of data to be processed, how the calculations were implemented, the compiler and the packages used, as well as the capabilities of the hardware employed.Thus, processing times may vary as any of these factors change.
The main advantages of this implementation are (1) the approach to obtain inversion solutions (parametric functional) is versatile, applicable to a wide variety of study objects because the linear regularizers that it integrates seek to promote general characteristics in the solutions, (2) the structural similarity regularizer allows reducing the uncertainty of the results, integrating two geophysical methods within the process of obtaining the solution, overcoming one of the great disadvantages of any geophysical method (ambiguity), finally (3) saves computational costs by reducing processing times which improves the efficiency of the process to obtain results and makes better use of available memory allowing to solve joint inversion problems for larger domains with less powerful hardware.
The use of other linear regularizers, such as the case of L, M, and N, can be incorporated without altering the performance too much since they are part of the expression I A (Equation ( 8)), which is calculated only once at the beginning of the inversion process.Therefore, changing the behavior of the solutions with other linear regularizers has no repercussions on the overall performance.
On the other hand, as can be seen in the convergence graphs (Figures 6, 10 and 13), the iterative process is unstable; that is, the models change with respect to the previous iteration by a high percentage.We associate this characteristic to the method we employ to minimize the parametric functional but overall to the use of many regularizers within the minimization process, as we noticed solutions with fewer regularizers active tend to be more stable, however, as also shown in the same figures, there is a tendency to stabilize the solution as the number of iterations increases, while in the RMS graph, it can be seen that the fit of the solutions is very good since early iterations.Therefore, this instability is not a problem for the quality of the results.
Finding the optimal combination of regularization parameter values where a solution exists and a low misfit of solutions with respect to observed data is achieved is a task that becomes more complex as more regularizers are non-zero in the parametric function.It is possible to resort to several methods, such as L-curve plots, which allow the moment of change between under-regularization and over-regularization to be visualized for many trials.The value of each regularization parameter controls the trade-off between fidelity to the observed data and the regularization conditions; in addition, other factors such as initial models, the quantity of a priori information, or data noise levels impact the inversion outcomes, making each particular case different.In this case, we overcame this drawback by using the trial-and-error approach, using near-zero values for regularization parameters and slowly increasing them until we achieved acceptable levels of error.Hence, several trials were performed to obtain these results (Figures 5 and 9).This is an opportunity area to improve by finding a way to obtain good results semi-automatically.

Figure 1 .
Figure 1.3D mesh illustrating the cells involved in the computation of discrete numerical derivatives depending on the position of the central node.

Figure 1 .
Figure 1.3D mesh illustrating the cells involved in the computation of discrete numerical derivatives depending on the position of the central node.

Figure 2 .
Figure 2. Radar chart ofTable1 normalized values, highlighting the advantages and disadvantages of each routine as the inversion domain increases.

Figure 3 .
Figure 3. Synthetic model of an inclined dike of 1 g/cm 3 in density and 1 A/m in magnetization buried at a depth of 50 m in a homogeneous bottom of zero contrast, (a) profile, (b) plan section at 50 m depth, (c) plan section at 400 m depth and (d) 3-D view of the buried dike.

Figure 3 .
Figure 3. Synthetic model of an inclined dike of 1 g/cm 3 in density and 1 A/m in magnetization buried at a depth of 50 m in a homogeneous bottom of zero contrast, (a) profile, (b) plan section at 50 m depth, (c) plan section at 400 m depth and (d) 3-D view of the buried dike.

Figure 3 .
Figure 3. Synthetic model of an inclined dike of 1 g/cm 3 in density and 1 A/m in magnetizati buried at a depth of 50 m in a homogeneous bottom of zero contrast, (a) profile, (b) plan section 50 m depth, (c) plan section at 400 m depth and (d) 3-D view of the buried dike.

Figure 4 .
Figure 4. (a) Gravity anomaly (mGal) and (b) magnetic anomaly (nT) produced by buried di model in Figure 3.The data were contaminated with Gaussian noise of 2% of the standard deviati in order to test the noise sensitivity of the inversion solutions.

Figure 4 .
Figure 4. (a) Gravity anomaly (mGal) and (b) magnetic anomaly (nT) produced by buried dike model in Figure 3.The data were contaminated with Gaussian noise of 2% of the standard deviation in order to test the noise sensitivity of the inversion solutions.
z cell and 1 × 10 8 rest of domain 0.3 top z cell and 1 × 10 8 rest of domain Std_Dev (2) m 0.3 top z cell and 1 × 10 8 rest of domain 0.3 top z cell and 1 × 10 8 rest of domain Std_Dev

Figure 5 .
Figure 5. Profiles of inversion results sliced at 500 m in  axis.Density distribution models (g/cm 3 ) for (a) standard separate inversion, (b) our robust joint inversion.Magnetization distribution (A/m) for (c) standard separate inversion and (d) our robust joint inversion.Black markers represent the original position of buried dike and red markers represent Well location used as a priori information.

Figure 6 .
Figure 6.Model convergence (Equation (18)) and data misfit (Equation (19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

Figure 5 .
Figure 5. Profiles of inversion results sliced at 500 m in y axis.Density distribution models (g/cm 3 ) for (a) standard separate inversion, (b) our robust joint inversion.Magnetization distribution (A/m) for (c) standard separate inversion and (d) our robust joint inversion.Black markers represent the original position of buried dike and red markers represent Well location used as a priori information.

Figure 5 .
Figure 5. Profiles of inversion results sliced at 500 m in  axis.Density distribution models (g/cm 3 ) for (a) standard separate inversion, (b) our robust joint inversion.Magnetization distribution (A/m) for (c) standard separate inversion and (d) our robust joint inversion.Black markers represent the original position of buried dike and red markers represent Well location used as a priori information.

Figure 6 .
Figure 6.Model convergence (Equation (18)) and data misfit (Equation (19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

Figure 6 .
Figure 6.Model convergence (Equation (18)) and data misfit (Equation (19)) were achieved in inversion of test data 1 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.

Figure 7 .
Figure 7. Synthetic test model 2 of simple geometric buried bodies with +1 g/cm 3 density contrast, +1 A/m magnetization contrast (rose bodies) and −1 g/cm 3 density contrast, −1 A/m magnetization contrast (blue body).(a) profile view at y = 250 m, (b) profile view at y = 700 m, (c) top view at z = 0 m and (d) 3D view.

Figure 8 .
Figure 8.(a) Gravity anomaly (mGal) and (b) magnetic anomaly (nT) produced by test model 2 in Figure 7.The data were contaminated with Gaussian noise of 5% of the standard deviation in order to test the noise sensitivity of the inversion solutions.

Figure 7 .
Figure 7. Synthetic test model 2 of simple geometric buried bodies with +1 g/cm 3 density contrast, +1 A/m magnetization contrast (rose bodies) and −1 g/cm 3 density contrast, −1 A/m magnetization contrast (blue body).(a) profile view at y = 250 m, (b) profile view at y = 700 m, (c) top view at z = 0 m and (d) 3D view.

Figure 7 .
Figure 7. Synthetic test model 2 of simple geometric buried bodies with +1 g/cm 3 density contrast, +1 A/m magnetization contrast (rose bodies) and −1 g/cm 3 density contrast, −1 A/m magnetization contrast (blue body).(a) profile view at y = 250 m, (b) profile view at y = 700 m, (c) top view at z = 0 m and (d) 3D view.

Figure 8 .
Figure 8.(a) Gravity anomaly (mGal) and (b) magnetic anomaly (nT) produced by test model 2 in Figure 7.The data were contaminated with Gaussian noise of 5% of the standard deviation in order to test the noise sensitivity of the inversion solutions.

Figure 8 .
Figure 8.(a) Gravity anomaly (mGal) and (b) magnetic anomaly (nT) produced by test model 2 in Figure 7.The data were contaminated with Gaussian noise of 5% of the standard deviation in order to test the noise sensitivity of the inversion solutions.

1 27 Figure 9 .
Figure 9. Profile comparation of separate vs. joint inversion results of test data 2. Density distribution models (g/cm 3 ) and Magnetization distribution (A/m) sliced at (a) 300 m in -axis and (b) 750 m in -axis to better show position of four buried bodies.Rose and blue markers represent the original position of buried bodies.

Figure 9 .
Figure 9. Profile comparation of separate vs. joint inversion results of test data 2. Density distribution models (g/cm 3 ) and Magnetization distribution (A/m) sliced at (a) 300 m in y-axis and (b) 750 m in y-axis to better show position of four buried bodies.Rose and blue markers represent the original position of buried bodies.

Figure 9 .
Figure 9. Profile comparation of separate vs. joint inversion results of test data 2. Density distribution models (g/cm 3 ) and Magnetization distribution (A/m) sliced at (a) 300 m in -axis and (b) 750 m in -axis to better show position of four buried bodies.Rose and blue markers represent the original position of buried bodies.

Figure 10 .
Figure 10.Model convergence (Equation (18)) and data misfit (Equation (19)) achieved in inversion of test data 2 at each iteration, (a) for standard separate inversion and (b) for our robust joint inversion.
2 nT absolute accuracy.The magnetic field in the area is characterized by B = 41,922.8,I = 50.86,and D = 4.51.Appl.Sci.2023, 13, x FOR PEER REVIEW 16 of 27 2 nT absolute accuracy.The magnetic field in the area is characterized by B = 41,922.8,I = 50.86,and D = 4.51.

Figure 11 .
Figure 11.Topography of the maar Los Contreras, SLP.Acquired gravity and magnetic data location in black marks.

Figure 11 .
Figure 11.Topography of the maar Los Contreras, SLP.Acquired gravity and magnetic data location in black marks.

Figure 12 .
Figure 12.Observed geophysical data in Los Contreras area.(a) Residual complete Bouguer gravity anomaly and (b) Residual magnetic anomaly.

Figure 12 .
Figure 12.Observed geophysical data in Los Contreras area.(a) Residual complete Bouguer gravity anomaly and (b) Residual magnetic anomaly.

27 Figure 13 .
Figure 13.Inversion results of preliminary field data, slicing the 3-D view each 300 meters from top to bottom until no more physical property contrast is visible, (a) Density contrast, (b) Magnetization contrast and (c) Convergence and RMS plots.

Figure 13 .
Figure 13.Inversion results of preliminary field data, slicing the 3-D view each 300 meters from top to bottom until no more physical property contrast is visible, (a) Density contrast, (b) Magnetization contrast and (c) Convergence and RMS plots.
1 × times faster than CPU (gfortran) version.Appl.Sci.2023, 13, x FOR PEER REVIEW 9 of 27 Figure 2. Radar chart of Table 1 normalized values, highlighting the advantages and disadvantages of each routine as the inversion domain increases.

Table 2 .
Inversion parameters for separate and joint inversion.

Table 3 .
Inversion parameters for separate and joint inversion.