# A prosperous and thorough analysis of gravity profiles for resources exploration utilizing the metaheuristic Bat Algorithm

Khalid S Essa, Omar A Gomaa, Mahmoud Elhussein, Yves Géraud, Marc Diraison, Zein E Diab

► To cite this version:

Khalid S Essa, Omar A Gomaa, Mahmoud Elhussein, Yves Géraud, Marc Diraison, et al.. A prosperous and thorough analysis of gravity profiles for resources exploration utilizing the metaheuristic Bat Algorithm. Scientific Reports, 2025, 15, pp.Article number: 5000. 10.1038/s41598-025-88350-4. hal-05033942

HAL Id: hal-05033942

https://hal.univ-lorraine.fr/hal-05033942v1

Submitted on 14 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A prosperous and thorough analysis of gravity profiles for resources exploration utilizing the metaheuristic Bat Algorithm

Khalid S. Essa<fnref n="1" />,²✉, Omar A. Gomaa<fnref n="1" />, Mahmoud Elhussein<fnref n="1" />, Yves Géraud², Marc Diraison² & Zein E. Diab<fnref n="1" />

Here, we present a remarkable methodology for unveiling subsurface structures with the potential to transform the exploration of mineral and ores resources, as well as the study of volcanic activity. By incorporating the Metaheuristic Bat algorithm (MBA) with the second horizontal gravity gradient (SHG) and employing variable window lengths, we aim to eliminate the regional effect in gravity data, thereby improving the precision of subsurface structure parameter estimation. Through rigorous evaluation on synthetic cases, we have demonstrated the robustness of our approach and its ability to handle diverse geological complexities and noise levels. Furthermore, our method has been applied to actual gravity data from three distinct locations: Canada, India, and Cuba, yielding excellent results that confirm the reliability and applicability of our methodology to real-world geological settings. We are confident that the use of variable window lengths in the SHG computation, coupled with the optimization of the global optimal solution via the Metaheuristic Bat Algorithm, can significantly contribute to the enhanced precision of subsurface structural parameter estimation. We hope our research will inspire others to explore this groundbreaking methodology and continue advancing the field of subsurface structure optimization.

**Keywords** MBA, Gravity data analysis, SHG, Source parameters estimation, Volcanic and mining studies

Gravity data serves as a powerful, non-invasive tool for exploring the Earth's subsurface, revealing the distribution of various rock types, density variations, and hidden structures such as faults, cavities, buried channels, and changes in lithology<fnref n="1" />⁻⁵. Gravity surveys have a wide range of applications across multiple fields. They are instrumental in identifying potential exploration targets and locating mineral and ore deposits. Additionally, these surveys are essential for studying the structure of volcanoes and magma chambers, which helps in assessing volcanic hazards⁶. In the context of climate change, gravity surveys are used to measure the thickness of glaciers and ice sheets, providing critical data for understanding environmental shifts. Furthermore, they determine the Earth's shape and its gravitational field, essential for accurate positioning systems like GPS. Finally, gravity surveys assist in locating vital groundwater resources, especially in arid regions⁷⁻¹³.

Unveiling the Earth's subsurface through gravity data analysis presents a formidable scientific challenge. Diverse methods have been developed to interpret this data and accurately image hidden subsurface structures. These methods can be broadly categorized as follows; Traditional approaches utilize well-established mathematical techniques to analyze gravity data. Forward modeling involves creating a theoretical subsurface model and calculating its gravitational response to compare with measured data¹⁴. Inversion techniques directly infer subsurface properties from gravity data, posing significant mathematical challenges¹⁵. Additionally, these methods often approximate the subsurface with simplified geometric shapes, such as spheres, cylinders, or faults, to aid in analysis¹⁶. However, more complex models incorporating intricate geological features require additional data and computational power, providing a more accurate representation of subsurface¹⁷.

Non-Traditional approaches explore advanced algorithms and techniques to extract meaningful insight from gravity data. These methods often integrate artificial intelligence, machine learning, and bio-inspired algorithms such as the Bat Algorithm¹⁸. By addressing complex optimization challenges that the conventional techniques may find difficult, these approaches significantly enhance the interpretive power of gravity data analysis. Over the past decade, metaheuristic algorithms inspired by natural process have gained substantial attention. Among

<footnote n="1">¹Geophysics Department, Faculty of Science, Cairo University, Giza 12613, Egypt. ²GeoRessources Laboratory, University of Lorraine, Nancy 54500, France. ✉email: khalid_sa_essa@cu.edu.eg</footnote>
these, the Bat algorithm, which mimics bats' echolocation strategies, has proven particularly effective in solving intricate optimization problems¹⁹,²⁰.

Optimization algorithms play crucial role in gravity data analysis by estimating the parameters of subsurface structures. These algorithms are generally classified into two categories: traditional and non-traditional (metaheuristic) methods. Traditional optimization techniques, such as Gradient Descent²¹, Newton-Raphson Method²², Simulated Annealing²³, and Conjugate Gradient Method²⁴, rely on mathematical rigor and gradient-based computations to locate optimal solutions. While they perform efficiently for straightforward and well-behaved problems, their effectiveness diminishes in the face of complex, multi-model landscapes typical of geophysical data. In contrast, non-traditional metaheuristic algorithms are inspired by natural and biological processes, offering significant advantages in exploring complex, high-dimensional search spaces. These algorithms, include the Hunger Games Search Algorithm²⁵,²⁶, Genetic Algorithms²⁷, and Particle Swarm Optimization²⁸, excel in addressing the challenges posed by irregular and intricate optimization problems encountered in gravity data analysis.

The Bat algorithm belongs to the class of non-traditional metaheuristic algorithms and is recognized for its dynamic and flexible nature in addressing complex optimization problems. Inspired by the echolocation behavior of bats, it introduces a novel strategy for global optimization and has demonstrated significant effectiveness in parameters estimation for subsurface structures¹⁹. Unlike traditional methods, the Bat algorithm offers distinct advantages, particularly its ability to balance exploration and exploitation phases efficiently. By dynamically adjusting its strategies during the early stages of optimization, it achieves faster convergence rate²⁹. A key strength of the Bat Algorithm lies in its dual capability as both a global and local optimizer, making it well-suited for handling multi-modal problems commonly encountered in geophysical studies³⁰. Its adaptability and robustness across diverse geological environments enhance its reliability, especially when dealing with noisy and complex datasets. Despite its advantages, the Bat Algorithm is not without limitations. One notable drawback is its computational intensity, particularly when applied to large-scale datasets or highly intricate geological scenarios. Additionally, the algorithm's performance is sensitive to parameter selection, necessitating careful tuning to achieve optimal results.

This paper presents a novel enhanced metaheuristic algorithm designed to optimize subsurface structure parameters using gravity data. A key innovative of this approach is the incorporation of second horizontal gradients (SHG), which effectively reduce the influence of regional data up to the first order, thereby improving the accuracy of the results. By integrating SHG with the metaheuristic Bat algorithm (MBA), we achieve precise estimations of subsurface structures. The subsurface features are modeled as simple geometric shapes, including spheres, horizontal cylinders, and vertical cylinders. This modeling strategy further minimizes the impact of regional data and enhances the precision of the outcomes. The proposed algorithm offers a transformative approach to subsurface exploration, promising significant advancements in the accurate delineation of subsurface features and contributing to the efficient development of subsurface resources.

The subsequent sections of this paper provide an in-depth exploration of our proposed methodology. The “Introduction” section offers a comprehensive overview of the significance of gravity data, with particular attention to the detailed insights provided by horizontal gradients. The “Methodology” section outlines the integration of gravity data, the second horizontal gradients, and the Bat algorithm to optimize subsurface structure parameters effectively. In the “Uncertainty Analysis” section, we evaluate the robustness and accuracy of the algorithm through validation against various theoretical models subjected to different noise levels. Finally, the “Results and Discussion” section highlights the algorithm’s performance using both synthetic and real-world field datasets, demonstrating its effectiveness and applicability in real-world scenarios.

## The methodology

To achieve accurate results in interpreting gravity data for subsurface model parameters, it is crucial to employ a robust inversion algorithm with comprehensive capabilities. This ensures precise assessments of key parameters such as the depth, location, and shape of subsurface structures.

### Forward modeling

The gravity measurement at a specific location (x_j) along a survey profile (Fig. 1) can be expressed as³¹⁻³³,

$$g_{\text{total}}(x_j) = g_{\text{res}}(x_j, z, m, x_o, A, q) + g_{\text{reg}}(x_j, x_o) \quad (1)$$

$$g_{\text{res}}(x_j, z, m, x_o, A, q) = \frac{A^*z^m}{[(x_j - x_o)^2 + z^2]^q}, \quad j = 1, 2, 3, \dots, n \quad (2)$$

$$g_{\text{reg}}(x_j, x_o) = a(x_j - x_o) + b \quad (3)$$

where $g_{\text{total}}(x_j)$ is the total measured gravity anomaly, $g_{\text{res}}(x_j, z, m, x_o, A, q)$ is the residual anomaly, $g_{\text{reg}}(x_j, x_o)$ is the regional anomaly, (a, b) are constant values, $x_j$ and $x_o$ signify the measurement points, in meters, and the origin location of the target, z represents the depth, in meters, of the buried source, q is a dimensionless shape factor that describes the geometry of the buried mass, and A is the amplitude coefficient (mGal.m²q⁻¹), which depend on both shape (q) and another factor (m), yielding g in milligals (mGal), γ is the universal gravitational constant is (6.67384*10⁻¹¹) m³.kg⁻¹.s⁻², σ represents the density contrast (gm/cc), and r represents the body radius (m). Table 1 presents the different cases of subsurface bodies and their properties.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of three distinct geometric shape models. Panel (a) depicts a model of a vertically oriented cylinder model, Panel (b) presents a spherical model, and Panel (c) exhibits a horizontally positioned cylinder model.</figcaption></figure>

<table><caption>Table 1. Definition of A, q, and μ parameters for the different simple model cases.</caption><thead><tr><td>Case</td><td>μ</td><td>A</td><td>q</td></tr></thead><tbody><tr><td>Horizontal cylinder</td><td>1</td><td>2πγσ r<sup>2</sup></td><td>1</td></tr><tr><td>Vertical cylinder</td><td>0</td><td>πγσ r<sup>2</sup></td><td>0.5</td></tr><tr><td>Sphere</td><td>1</td><td>(4/3) πγσ r<sup>3</sup></td><td>1.5</td></tr></tbody></table>

## Second horizontal gradient (SHG)

The second horizontal gradient (SHG) technique is an essential tool employed to enhance gravity data for the estimation of subsurface structure parameters. This method is particularly important for highlighting nuanced variations in the gravity field, a concept widely discussed in geophysical literature, which underscores its role in improving the resolution of subsurface features³⁴. The use of variable window lengths in computing the SHG allows for adaptation to diverse geological conditions, which significantly improves the quality of data used for subsurface structures estimation. This approach effectively mitigates the influence of regional data and sharpens the resolution of nuanced geological attributes. Asfahani and Tlas³⁵ have corroborated this approach in their studies, demonstrating its efficacy in improving data quality.

To remove the regional background, the SHG operator was applied to Eq. (1). For three observation points along the gravity profile (x<sub>j</sub> - 2s, x<sub>j</sub>, x<sub>j</sub> + 2s), the SHG (g<sub>xx</sub> (x<sub>j</sub>, s)) can be expressed as described by Essa and Elhussein³⁶:

$$g_{xx}(x_j, s) = \frac{g_{\text{total}}(x_j + 2s) - 2g_{\text{total}}(x_j) + g_{\text{total}}(x_j - 2s)}{4s^2} \quad (4)$$
$$g_{xx}(x_j, s) = \left( \frac{\left( (A^*z^m / [(x_j - x_o + 2s)^2 + z^2]^q) + (a(x_j - x_o + 2s) + b) \right) \dots - 2 \left( (A^*z^m / [(x_j - x_o)^2 + z^2]^q) + (a(x_j - x_o) + b) \right) \dots}{4s^2} \right) + \left( (A^*z^m / [(x_j - x_o - 2s)^2 + z^2]^q) + (a(x_j - x_o - 2s) + b) \right) \dots \quad (5)$$

where s = 1, 2, 3, ..., N separation units are window lengths and x_j is the observation data point.

The optimized data acquired through the SHG, employing variable window lengths, serves as input for the Bat algorithm's optimization process. This integration is crucial for refining the estimation of subsurface structure parameters.

## Metaheuristic bat algorithm (MBA)

Yang³⁷ introduced the Metaheuristic Bat algorithm (MBA), drawing inspiration from the echolocation actions of micro-bats. These bats use echolocation to navigate and hunt in the dark, emitting noisy sound pulses within a range of 8 to 10 kHz and listening to the echoes bouncing off objects nearby. The Bat algorithm utilizes this echolocation behavior as a basis for optimizing objective functions.

The Bat algorithm operates through three key stages, each critical for effective optimization. In the first stage, bats use echolocation to measure distances, establishing a baseline for navigating the search space. In the second stage, bats fly at a consistent frequency within a specified range [Q_min, Q_max], starting with an initial velocity (V_i) and position (X_i), to locate target objects. This stage mimics the exploratory phase where bats search for optimal solutions. In the third stage, the loudness (L_i) and pulse emission rate (r_i) dynamically adjust based on their proximity to the target, enhancing the exploitation phase where fine-tuning occurs. The frequency range [Q_min, Q_max] corresponds to the wavelength spectrum [K_min, K_max], which can be modified to alter the bats' movement range, optimizing their search behavior as per Eqs. (6-8). The algorithm updates loudness and emission rates only if new solutions indicate improvement, signaling progress towards the optimal solution³⁷ Eqs. (9, 10). This study investigates the impact of optimizing parameters such as frequency (Q_i), loudness (L_i), and pulse rate (r_i) on the Bat Algorithm's convergence rate, exploring various parameter ranges to identify optimal settings for enhanced performance.

Figure 2 illustrates the impact of each tuning parameter set (Q_i, L_i, and r_i) on convergence behavior, highlighting the optimal set as (Q_i = [0, 5], L_i = 1.0, and r_i = 0.9). This configuration demonstrates the lowest Normalized Root Mean Square Error (NRMSE) for the objective function and achieves rapid convergence to the optimal solution. The pulse rate (r_i) ranges from 0 to 1, representing no pulse emission to maximum emission,

<figure><img src="image_4.png" /><figcaption>Fig. 2. The convergence curves associated with the solution of optimal parameters.</figcaption></figure>
depending on target proximity. Initial loudness ($L_i$) typically ranges between 1 and 2, decreasing as bats approach prey, while pulse emission rate increases. Initially, the initial speed ($V_i$) at position ($X_i$) is set to zero at the beginning of the Bat algorithm inversion process. Selecting an appropriate frequency or wavelength range is crucial, ideally aligning with the domain of awareness before narrowing it down. In this study, a frequency range of [0, 5] was identified as optimal after testing various values.

The mention-below equations described the relationship between algorithm parameters, as outlined by Yang³⁷.

$$Q_i^t = Q_{min} + (Q_{max} - Q_{min}) * \beta \qquad (6)$$

$$V_i^{(t+1)} = V_i^t + (X_i^t - X_{best}) * Q_i^t \qquad (7)$$

$$X_i^{(t+1)} = X_i^t + V_i^{(t+1)} \qquad (8)$$

$$L_i^{(t+1)} = \alpha L_i^t \qquad (9)$$

$$r_i^t = r_i^o [1 - \exp(-\gamma \tau)] \qquad (10)$$

where, $Q_i$ represents the frequency of the ith bat, which is updated in each iteration, $\beta$ signifies a random vector sample from uniform distribution within the range [0, 1], and $X_{best}$ denotes currently the global best solution among all bat numbers, $\alpha$ and $\gamma$ are constants, where $0 < \alpha < 1$ and $\gamma > 0$ and $\tau$ is the scaling factor.

The Bat algorithm generates new results from each chosen best solution in the local search using a random path, as shown below:

$$X_{new} = X_{old} + \epsilon L^t \qquad (11)$$

where $\epsilon$ is a random number within [-1, 1], and $L^t$ characterizes the average loudness of all bat numbers at the existing stage, $t$ represents the number of iterations. In terms of accuracy and performance, the MBA outperforms most other algorithms. Essentially, if the frequency variations are replaced by a random parameter and $L_i=1$ and $r_i=1$ are set, the bat algorithm effectively becomes the regular particle swarm optimization (PSO)³⁷.

A new code has been developed called the “Inversion MBA”, which offers a promising solution for identifying the optimal subsurface sources that accurately represent the real data. This code focuses on significant parameters such as depth, location, body shape, and amplitude coefficient ($z, x_o, q, A, m$) to ensure an optimal solution is reached. The iterative process involves bats moving randomly within the search space, refining their solutions and determining the optimal positions based on solution evaluations. The best solution, characterized by the lowest misfit function value, is denoted as the best location ($X_{best}$). In each iteration, the best solutions are compared, and the optimal solution is selected. The final best solution ($X_{best}$) is picked after a specified number of iterations. The MBA inversion code has undergone thorough testing on various synthetic models before being applied to real datasets, ensuring accurate and reliable results.

In our novel approach to subsurface parameter estimation, we introduce an algorithm that seamlessly integrates second horizontal gradients (SHG) using varying window values. The process unfolds as follows: We initiate by generating objective functions based on gravity data sets. These functions serve as fitness criteria for evaluating the fit between observed and calculated gravity anomalies. Next, we meticulously apply SHG within distinct windows, ensuring robustness across anomalies. The SHG, represented by Eq. (4), enhances the sensitivity of anomaly detection, capturing nuanced variations in the data. Each SHG anomaly becomes a potential target for optimization. Our secret weapon lies in the Bat algorithm (BA), inspired by the echolocation behavior of bats. The flowchart detailing this algorithm is visually represented in Fig. 3. Here’s how it works: We initialize virtual bats within the search space. Each bat corresponds to a potential solution. The parameters of the bat algorithm include position ($X_i$), variables standing for subsurface parameters (e.g., depth ($z$), origin ($x_o$), amplitude coefficient ($A$), shape factors ($q$), and ($m$)). These are randomly selected from the search space. Frequencies ($Q_i$), guiding the bat’s movement during optimization. Velocities ($V_i$), controlling exploration speed. Loudness ($L_i$), influencing the bat’s behavior. Pulse rates ($r_i$), introducing randomness.

The objective function ($F_{obj}$), defined by Eq. (12), and quantifies the fit. It calculates the normal root-mean-square error (NRMSE) between observed and calculated gravity anomalies. The bat with the least misfit becomes ($X_{best}$). While the program hasn’t reached the maximum iterations: Modify frequencies (Eq. 6), and update velocities and positions (Eqs. 7, 8) to create fresh solutions, if a random number exceeds ($r_i$) the algorithm selects a solution from the best set. Then Generates a local solution around it (Eq. 11). Accept the fresh solutions after comparing the current solution with the initial ones, adjust ($r_i$) and ($L_i$) (Eqs. 9, 10), rank the bats that stand for the solutions and decide the best solution ($X_{best}$).

After applying BA to different gradient anomalies, we aggregate the results. Calculating their average significantly reduces uncertainty and minimizes parameter estimation errors. By leveraging the collective intelligence of BA and statistical averaging, our algorithm enhances accuracy in identifying subsurface features. The resulting parameter estimates exhibit improved reliability, making them invaluable for geophysical exploration.

$$F_{obj} = \frac{100}{\text{Max}(g_{\text{obs}}) - \text{Min}(g_{\text{obs}})} \sqrt{\frac{\sum_{j=1}^{N} [g_{\text{obs}} - g_{\text{cal}}]^2}{N}} \qquad (12)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. The flowchart outlines the procedure for the current MBA method.</figcaption></figure>
where N is the data points numbers, $g_{\text{obs}}$ represent the observed gravity and $g_{\text{cal}}$ signifies the calculated gravity. Initially, Eq. (12) is utilized to calculate misfits, and then the bat with the least misfit is picked as the $X_{\text{best}}$.

## Results and discussion
### Synthetic datasets

The effectiveness of the proposed MBA approach for inverting gravity data under various challenging conditions was thoroughly evaluated. We created synthetic datasets with different complexities, including varying noise levels, background (regional) effects, and neighboring effects.

By analyzing and comparing the MBA's performance against established models on these complex datasets, we conducted a comprehensive assessment of its ability to manage data imperfections. The robustness and reliability of the MBA were rigorously tested, both in terms of handling the original data and in the accuracy of the predicted datasets it generates.

### Model 1: effect of various noise levels

This study investigates the effectiveness of the MBA scheme in recovering source parameters from simulated gravity anomaly data.

We first generated a theoretical gravity anomaly using a specific set of parameters ($A=140$ mGal.km, $z=7$ km, $q=1.0$, $x_0 = 0$ km, and $m=1$). This anomaly was calculated over a 101-kilometer profile length using Eq. (2) (Fig. 4a). The resulting data served as a baseline for the analysis. The MBA technique was then applied to optimize the inversion process and obtain the most accurate representation of the gravity response. This involved analyzing second horizontal gradient anomalies (SHG) for various values of parameter 's' ($s=1, 2, 3, 4, 5, 6, 7, 8, 9, 10$ km) (Fig. 4b). The optimization aimed to minimize the normalized root mean square error (NRMSE) of the objective function. Figure 4c depicts the evolution of the minimum NRMSE, signifying the best solution overall, across the total number of iterations. The plot shows that the NRMSE reaches its minimum value after approximately 50 iterations for all simulated "microbats" (individual elements within the MBA). Additionally, Fig. 4d illustrates the average NRMSE for each iteration, indicating convergence towards the optimal solution. The successful reduction in the NRMSE suggests the efficacy of the MBA method in retrieving the original parameters. This is supported by the evaluation of the generated gravity data. Table 2 displays the average

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** First model, a horizontal cylinder, in case of noise-free conditions. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{\text{Obj}}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>
<table><caption>**Table 2.** Results of parameter estimation for the noise-free horizontal cylinder synthetic model.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A (mGal. km)</td><td>140</td><td>10-500</td><td>140</td><td>140</td><td>138</td><td>140</td><td>140</td><td>140</td><td>142</td><td>140</td><td>140</td><td>140</td><td>140 ± 0.94</td><td rowspan="4">0.000</td></tr><tr><td>z (km)</td><td>7</td><td>1-10</td><td>7</td><td>7</td><td>7</td><td>7.1</td><td>7</td><td>7</td><td>6.9</td><td>7</td><td>7</td><td>7</td><td>7 ± 0.047</td></tr><tr><td>x<sub>o</sub> (km)</td><td>0</td><td>-10-10</td><td>0</td><td>0</td><td>0</td><td>0.1</td><td>0</td><td>0</td><td>0</td><td>-0.1</td><td>0</td><td>0</td><td>0 ± 0.04</td></tr><tr><td>q</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td></tr><tr><td>m</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** First model, a horizontal cylinder, in case of a 5% noise level condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Ob}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

estimated parameters derived from various SHG anomalies. The parameters are as follows: $A = 140 \pm 0.94$ mGal. km, $z = 7 \pm 0.047$ km, $x_o = 0 \pm 0.04$ km, $q = 1 \pm 0$, $m = 1 \pm 0$, these parameters collectively characterize a horizontal cylinder. The normalized root mean square error (NRMSE) was found to be 0.000.

To assess the robustness of the MBA technique, we introduced noise into the synthetic data at two levels: 5% and 15%.

At a 5% noise level (depicted in Fig. 5), the MBA successfully identified the optimal model parameters by minimizing the NRMSE of the goal function (Fig. 5a). The estimated parameters closely resembled the original values. Figure 5b displays the SHG anomalies for this case. Figure 5c and d visually represent the convergence process through minimum and mean NRMSE values.

Subsequently, Table 3 outlines the results of the estimated parameters post our algorithm's application to the model with a 5% noise level. The estimated parameters are as follows: $A = 145 \pm 4.55$ mGal.km, $z = 6.9 \pm 0.14$ km, $x_o = 0 \pm 0.35$ km, $q = 1 \pm 0$, $m = 1 \pm 0$. The NRMSE was determined to be 0.00015.
<table><caption>**Table 3.** Results of parameter estimation for the horizontal cylinder synthetic model with a 5% noise level.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A (mGal. km)</td><td>140</td><td>10-500</td><td>138</td><td>140</td><td>142</td><td>144</td><td>148</td><td>146</td><td>148</td><td>150</td><td>152</td><td>142</td><td>145 ± 4.55</td><td rowspan="4">0.00015</td></tr><tr><td>z (km)</td><td>7</td><td>1-10</td><td>6.8</td><td>7</td><td>6.9</td><td>7</td><td>7</td><td>6.8</td><td>7.1</td><td>7</td><td>6.7</td><td>6.7</td><td>6.9 ± 0.14</td></tr><tr><td>x<sub>o</sub> (km)</td><td>0</td><td>-10-10</td><td>0.2</td><td>0.6</td><td>0</td><td>-0.2</td><td>0.4</td><td>-0.4</td><td>0</td><td>0</td><td>-0.6</td><td>0</td><td>0 ± 0.35</td></tr><tr><td>q</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td></tr><tr><td>m</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td><td></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** First model, a horizontal cylinder, in case of a 15% noise level condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution (F<sub>Obj</sub>) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

Similarly, at a 15% noise level (depicted in Fig. 6), the MBA approach yielded reasonable results (Fig. 6a). The estimated parameters again showed good agreement with the originals. The SHG anomalies are displayed in Fig. 6b, and convergence plots (Fig. 6c and d) support these findings. Table 4 presents the results of the estimated parameters after applying our algorithm to the model with a 15% noise level. The parameters are: A = 130 ± 5.77 mGal.km, z = 7.2 ± 0.23 km, x<sub>o</sub> = 1 ± 1.05 km, q = 1 ± 0, m = 1 ± 0. The NRMSE is 0.00019. Figure 7 shows the relative errors for each model parameter, showcasing the robustness and accuracy of our algorithm across under varying noise levels. The results highlight the algorithm's ability to maintain reliable performance, even in the presence of substantial noise. These findings emphasize the effectiveness of the MBA technique in accurately estimating subsurface parameters, further validating its reliability and robustness in noisy environments.

**Model 2: effect of regional anomaly**

This section evaluates the MBA performance under more challenging conditions. We introduced a more complex scenario by combining: A numerical model with specific parameters (A = 600 mGal.km², z = 2 km,
<table><caption>**Table 4.** Results of parameter estimation for the horizontal cylinder synthetic model with a 15% noise level.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s=1 km</th><th>s=2 km</th><th>s=3 km</th><th>s=4 km</th><th>s=5 km</th><th>s=6 km</th><th>s=7 km</th><th>s=8 km</th><th>s=9 km</th><th>s=10 km</th></tr></thead><tbody><tr><td>A (mGal. km)</td><td>140</td><td>10-500</td><td>130</td><td>125</td><td>135</td><td>140</td><td>130</td><td>120</td><td>125</td><td>130</td><td>130</td><td>135</td><td>130 ± 5.77</td><td rowspan="4">0.00019</td></tr><tr><td>z (km)</td><td>7</td><td>1-10</td><td>7.1</td><td>6.8</td><td>7.4</td><td>7.2</td><td>7.1</td><td>7.4</td><td>6.9</td><td>7.5</td><td>7.2</td><td>7.4</td><td>7.2 ± 0.23</td></tr><tr><td>x<sub>o</sub> (km)</td><td>0</td><td>-10-10</td><td>0</td><td>1</td><td>2</td><td>1</td><td>-1</td><td>2</td><td>0</td><td>2</td><td>1</td><td>2</td><td>1 ± 1.05</td></tr><tr><td>q</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td></tr><tr><td>m</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td><td></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** The robustness and accuracy of the parameter estimations for the first model under varying noise levels.</figcaption></figure>

$q=1.5$, $x_o = 0$ km, and $m=1$) collectively characterize a sphere body calculated using Eq. (1), a linear regional anomaly ($1.5*(x_j - x_o) + 25$), and 10% random noise added to the data (Figs. 8 and 9).

Despite the added complexity, the MBA technique was employed to estimate the original spherical parameters from the SHG anomalies for various 's' values ($s=1$ to $10$ km) (Fig. 8b). The results were encouraging. The MBA successfully recovered parameters close to the originals: $A = 600 \pm 9.43$ mGal.km$^2$, $z = 1.9 \pm 0.13$ km, $q = 1.5 \pm 0.05$, $x_o = 0 \pm 0$ km, and $m = 1 \pm 0$, with a normalized root mean square error (NRMSE) of 0.00027 (Table 5). This indicates a good match between the calculated and the anomaly (Fig. 8a). Similar to previous tests, Fig. 8c shows the optimization process. It depicts the evolution of the minimum NRMSE, representing the best overall solution, across the total number of iterations. As expected, the NRMSE gradually decreases, reaching its minimum after approximately 50 iterations for all “microbats” within the MBA. Figure 8d complements this by illustrating the average NRMSE for each iteration, further indicating convergence towards the optimal solution.

After introducing a 10% noise level, the MBA was once again utilized to derive the spherical parameters. The resulting estimates were as follows: $A = 595 \pm 30.46$ mGal.km$^2$, $z = 2 \pm 0.28$ km, $q = 1.4 \pm 0.13$, $x_o = 0 \pm 2.31$ km, and $m = 1 \pm 0.12$, accompanied with a NRMSE of 0.00014 (Table 6). This suggests a strong agreement between the computed and the noisy anomaly (Fig. 9a). The SHG anomalies are displayed in Fig. 9b. In line with the original investigation, Fig. 9c portrays the optimization process, illustrating the diminishing NRMSE over iterations and converging towards the optimal solution. As anticipated the NRMSE gradually decreases achieving its minimum after roughly 50 iterations across all MBA “microbats”. Figure 9d complements this by
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Second model, a sphere, with a first-order regional and in case of noise-free condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying ‘s’ window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

illustrating the average NRMSE for each iteration, further indicating convergence towards the optimal solution. These findings demonstrate the robustness of the MBA method, even when dealing with increased noise levels. Figure 10 presents the relative errors for each model parameter, demonstrating the algorithm's capability to accurately estimate subsurface parameters even in the presence of regional anomalies and noise. This highlights the robustness and precision of the proposed approach in challenging conditions.

## Model 3: Effect of neighboring structures

To evaluate the stability and efficacy of the proposed algorithm, two scenarios were studied. The first scenario involved two separated vertical cylinder models. The first model consisted of a body with $A_1 = 230$ mGal.km, $z_1 = 5$ km, $q_1 = 0.5$, $x_{o1} = -20$ km, and $m_1 = 0$, while the second model consisted of a body with $A_2 = 200$ mGal.km, $z_2 = 3$ km, $q_2 = 0.5$, $x_{o2} = 25$ km, and $m_2 = 0$, these parameters collectively characterize a vertical cylinder (Fig. 11a) were computed using Eq. (2) and contaminated with 10% random noise (Fig. 12a). The MBA method was used to estimate the model parameters for all SHG anomalies (Figs. 11b and 12b). Table 7, the results show that the estimated parameters for the first body are $A_1 = 225 \pm 10.80$ mGal.km, $z_1 = 5.2 \pm 0.30$ km, $q_1 = 0.5 \pm 0.03$, $x_{o1} = -20 \pm 1.33$ km, and $m_1 = 0 \pm 0$, the second body are $A_2 = 205 \pm 4.71$ mGal.km, $z_2 = 3 \pm 0.18$ km, $q_2 = 0.5 \pm 0.05$, $x_{o2} = 25 \pm 0.94$ km, and $m_2 = 0 \pm 0$, respectively. The NRMSE-misfit is 0.0000023, and the predicted and original anomalies comparison is displayed in Fig. 11a. Similarly, Fig. 11c shows the optimization process. The minimum NRMSE, representing the best solution overall, decreases steadily across iterations. This indicates convergence towards the optimal parameters, achieved after approximately 160 iterations for all “microbats” within the MBA. Figure 11d complements this by illustrating the average NRMSE per iteration.

With the introduction of a 10% noise level, the estimated parameters for the first body ($A_1 = 240 \pm 18.26$ mGal.km, $z_1 = 4.8 \pm 0.44$ km, $q_1 = 0.5 \pm 0.06$, $x_{o1} = -19 \pm 2.11$ km, and $m_1 = 0 \pm 0$) and the second body ($A_2 = 200 \pm 6.67$ mGal.km, $z_2 = 3.4 \pm 0.15$ km, $q_2 = 0.5 \pm 0.06$, $x_{o2} = 25 \pm 1.05$ km, and $m_2 = 0 \pm 0$) are presented in Table 8. The resulting NRMSE-misfit of 0.0000032 indicates a good match between the predicted and original anomalies, despite the additional noise (Fig. 12a). The optimization process depicted in Fig. 12c, demonstrates a consistent decrease in the minimum NRMSE across iterations, signifying convergence towards the optimal parameters after approximately 260 iterations for all “microbats” within the MBA. Figure 12d supplements this by illustrating the average NRMSE per iteration, further confirming the convergence towards the optimal solution.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Second model, a sphere, with a first-order regional and in case of a 10% noise level condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

<table><caption>**Table 5.** Results of parameter estimation for the spherical synthetic model with a first order regional in the absence of noise.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A (mGal. km<sup>2</sup>)</td><td>600</td><td>300-800</td><td>600</td><td>600</td><td>590</td><td>610</td><td>590</td><td>600</td><td>620</td><td>600</td><td>590</td><td>600</td><td>600 ± 9.43</td><td></td></tr><tr><td>z (km)</td><td>2</td><td>1-10</td><td>1.7</td><td>1.8</td><td>2</td><td>1.9</td><td>2.1</td><td>2</td><td>2</td><td>1.8</td><td>1.9</td><td>1.8</td><td>1.9 ± 0.13</td><td rowspan="2">0.00027</td></tr><tr><td>x<sub>o</sub> (km)</td><td>0</td><td>-10-10</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0 ± 0</td></tr><tr><td>q</td><td>1.5</td><td>0-2</td><td>1.5</td><td>1.4</td><td>1.6</td><td>1.5</td><td>1.5</td><td>1.45</td><td>1.5</td><td>1.55</td><td>1.5</td><td>1.5</td><td>1.5 ± 0.05</td><td></td></tr><tr><td>m</td><td>1</td><td>0-2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1 ± 0</td><td></td></tr></tbody></table>

<table><caption>**Table 6.** Results of parameter estimation for the spherical synthetic model with a first order regional with a 10% noise level.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A (mGal. km<sup>2</sup>)</td><td>600</td><td>300-800</td><td>580</td><td>570</td><td>640</td><td>630</td><td>585</td><td>615</td><td>630</td><td>575</td><td>565</td><td>560</td><td>595 ± 30.46</td><td></td></tr><tr><td>z (km)</td><td>2</td><td>1-10</td><td>1.8</td><td>1.9</td><td>2.4</td><td>2.3</td><td>1.8</td><td>1.6</td><td>2</td><td>2.3</td><td>1.7</td><td>2.2</td><td>2 ± 0.28</td><td rowspan="2">0.00014</td></tr><tr><td>x<sub>o</sub> (km)</td><td>0</td><td>-10-10</td><td>-3</td><td>2</td><td>1</td><td>-3</td><td>0</td><td>4</td><td>-1</td><td>2</td><td>0</td><td>-2</td><td>0 ± 2.31</td></tr><tr><td>q</td><td>1.5</td><td>0-2</td><td>1.3</td><td>1.4</td><td>1.35</td><td>1.55</td><td>1.3</td><td>1.2</td><td>1.5</td><td>1.65</td><td>1.35</td><td>1.4</td><td>1.4 ± 0.13</td><td></td></tr><tr><td>m</td><td>1</td><td>0-2</td><td>0.8</td><td>1.2</td><td>1.1</td><td>1</td><td>0.9</td><td>1.1</td><td>1</td><td>0.9</td><td>1</td><td>1</td><td>1 ± 0.12</td><td></td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** The robustness and accuracy of the parameter estimations for the second model, in case of existence of regional anomaly without and with 10% random noise.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Third model, two separated structures, in case of noise-free conditions. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

In the second scenario, the structures were repositioned closer together. The first model parameters were $A_1 = 230$ mGal.km, $z_1 = 5$ km, $q_1 = 0.5$, $x_{o1} = -10$ km, and $m_1 = 0$, while the second model parameters were $A_2 = 200$ mGal.km, $z_2 = 3$ km, $q_2 = 0.5$, $x_{o2} = 10$ km, and $m_2 = 0$ (Fig. 13a). Computing using Eq. (2) and contaminated with 10% random noise (Fig. 14a). The MBA method was used to estimate the model parameters for all SHG anomalies (Figs. 13b and 14b). Table 9, the results show that the estimated parameters for the first body are $A_1 = 235 \pm 11.55$ mGal.km, $z_1 = 5.1 \pm 0.33$ km, $q_1 = 0.5 \pm 0.02$, $x_{o1} = -10 \pm 1.56$ km, and $m_1 = 0.011 \pm 0.01$, the second body are $A_2 = 195 \pm 8.82$ mGal.km, $z_2 = 2.9 \pm 0.18$ km, $q_2 = 0.5 \pm 0.03$, $x_{o2} = 10 \pm 1.16$ km, and $m_2 = 0.01 \pm 0.01$, respectively. The NRMSE-misfit is 0.00000215, and the predicted and original anomalies comparison is displayed in Fig. 13a. Similarly, Fig. 13b and c shows the SHG anomalies and the optimization process. The minimum
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Third model, two separated structures, in case of a 10% noise level condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

<table><caption>**Table 7.** Results of parameter estimation for the two separated synthetic model under a noise-free condition.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A<sub>1</sub> (mGal.km)</td><td>230</td><td>100-600</td><td>230</td><td>210</td><td>215</td><td>235</td><td>235</td><td>225</td><td>230</td><td>240</td><td>220</td><td>210</td><td>225 ± 10.80</td><td rowspan="10">0.0000023</td></tr><tr><td>z<sub>1</sub> (km)</td><td>5</td><td>1-10</td><td>5.4</td><td>5.6</td><td>4.8</td><td>5</td><td>5.2</td><td>5</td><td>4.8</td><td>5.2</td><td>5.4</td><td>5.6</td><td>5.2 ± 0.30</td></tr><tr><td>x<sub>o1</sub> (km)</td><td>-20</td><td>-30-30</td><td>-20</td><td>-18</td><td>-20</td><td>-20</td><td>-23</td><td>-19</td><td>-21</td><td>-19</td><td>-20</td><td>-20</td><td>-20 ± 1.33</td></tr><tr><td>q<sub>1</sub></td><td>0.5</td><td>0-2</td><td>0.5</td><td>0.55</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.55</td><td>0.45</td><td>0.5</td><td>0.5</td><td>0.45</td><td>0.5 ± 0.03</td></tr><tr><td>m<sub>1</sub></td><td>0</td><td>0-2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0 ± 0</td></tr><tr><td>A<sub>2</sub> (mGal.km)</td><td>200</td><td>100-600</td><td>205</td><td>200</td><td>210</td><td>200</td><td>200</td><td>205</td><td>205</td><td>215</td><td>205</td><td>205</td><td>205 ± 4.71</td></tr><tr><td>z<sub>2</sub> (km)</td><td>3</td><td>1-10</td><td>3</td><td>2.9</td><td>2.8</td><td>3</td><td>3.2</td><td>3</td><td>2.9</td><td>2.8</td><td>3.4</td><td>3</td><td>3 ± 0.18</td></tr><tr><td>x<sub>o2</sub> (km)</td><td>25</td><td>-30-30</td><td>24</td><td>25</td><td>24</td><td>26</td><td>24</td><td>25</td><td>25</td><td>26</td><td>24</td><td>27</td><td>25 ± 1.05</td></tr><tr><td>q<sub>2</sub></td><td>0.5</td><td>0-2</td><td>0.5</td><td>0.5</td><td>0.4</td><td>0.45</td><td>0.5</td><td>0.6</td><td>0.55</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5 ± 0.05</td></tr><tr><td>m<sub>2</sub></td><td>0</td><td>0-2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0 ± 0</td></tr></tbody></table>

NRMSE, representing the best solution overall, decreases steadily across iterations. This indicates convergence towards the optimal parameters, achieved after approximately 150 iterations for all “microbats” within the MBA. Figure 13d complements this by illustrating the average NRMSE per iteration.

With the introduction of a 10% noise level, the estimated parameters for the first body ($A_1 = 225 \pm 11.547$ mGal.km, $z_1 = 5.3 \pm 0.094$ km, $q_1 = 0.49 \pm 0.022$, $x_{o1} = -11 \pm 1.247$ km, and $m_1 = 0.01 \pm 0.008$) and the second body ($A_2 = 210 \pm 8.819$ mGal.km, $z_2 = 2.8 \pm 0.183$ km, $q_2 = 0.5 \pm 0.025$, $x_{o2} = 10 \pm 1.633$ km, and $m_2 = 0.005 \pm 0.005$) are presented in Table 10. The resulting NRMSE-misfit of 0.0000116 indicates a good match between the predicted and original anomalies, despite the additional noise (Fig. 14a). The optimization process depicted in Fig. 14c, demonstrates a consistent decrease in the minimum NRMSE across iterations, signifying convergence towards the optimal parameters after approximately 80 iterations for all “microbats” within the MBA. Figure 14d supplements this by illustrating the average NRMSE per iteration, further confirming the convergence towards the optimal solution. Figure 15 highlights the relative errors for each model parameter, affirming the stability
<table><caption>**Table 8.** Results of parameter estimation for the two separated synthetic model under a scenario with a 10% noise level.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s=1 km</th><th>s=2 km</th><th>s=3 km</th><th>s=4 km</th><th>s=5 km</th><th>s=6 km</th><th>s=7 km</th><th>s=8 km</th><th>s=9 km</th><th>s=10 km</th></tr></thead><tbody><tr><td>A<sub>1</sub> (mGal.km)</td><td>230</td><td>100-600</td><td>215</td><td>205</td><td>250</td><td>230</td><td>240</td><td>255</td><td>260</td><td>250</td><td>255</td><td>240</td><td>240 ± 18.26</td><td rowspan="10">0.0000032</td></tr><tr><td>z<sub>1</sub> (km)</td><td>5</td><td>1-10</td><td>4.6</td><td>4.8</td><td>4.5</td><td>5.4</td><td>5.5</td><td>5.2</td><td>4.4</td><td>4.6</td><td>4.8</td><td>4.2</td><td>4.8 ± 0.44</td></tr><tr><td>x<sub>o1</sub> (km)</td><td>-20</td><td>-30-30</td><td>-17</td><td>-16</td><td>-22</td><td>-17</td><td>-19</td><td>-18</td><td>-19</td><td>-21</td><td>-22</td><td>-19</td><td>-19 ± 2.11</td></tr><tr><td>q<sub>1</sub></td><td>0.5</td><td>0-2</td><td>0.45</td><td>0.45</td><td>0.55</td><td>0.55</td><td>0.6</td><td>0.5</td><td>0.4</td><td>0.55</td><td>0.5</td><td>0.45</td><td>0.5 ± 0.06</td></tr><tr><td>m<sub>1</sub></td><td>0</td><td>0-2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0 ± 0</td></tr><tr><td>A<sub>2</sub> (mGal.km)</td><td>200</td><td>100-600</td><td>195</td><td>190</td><td>210</td><td>205</td><td>205</td><td>200</td><td>205</td><td>200</td><td>200</td><td>190</td><td>200 ± 6.67</td></tr><tr><td>z<sub>2</sub> (km)</td><td>3</td><td>1-10</td><td>3.5</td><td>3.2</td><td>3.4</td><td>3.2</td><td>3.6</td><td>3.5</td><td>3.3</td><td>3.4</td><td>3.3</td><td>3.6</td><td>3.4 ± 0.15</td></tr><tr><td>x<sub>o2</sub> (km)</td><td>25</td><td>-30-30</td><td>25</td><td>24</td><td>24</td><td>25</td><td>25</td><td>25</td><td>26</td><td>24</td><td>25</td><td>27</td><td>25 ± 0.94</td></tr><tr><td>q<sub>2</sub></td><td>0.5</td><td>0-2</td><td>0.55</td><td>0.6</td><td>0.5</td><td>0.5</td><td>0.55</td><td>0.45</td><td>0.45</td><td>0.55</td><td>0.45</td><td>0.4</td><td>0.5 ± 0.06</td></tr><tr><td>m<sub>2</sub></td><td>0</td><td>0-2</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0 ± 0</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Third model, two neighboring structures, in case noise-free condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

and accuracy of the algorithm in handling complex geological conditions, even with closely spaced structures and noise.

## Uncertainty analysis

To evaluate the robustness and accuracy of our algorithm, extensive testing was conducted using three different theoretical models before applying it to real-world field datasets. This comprehensive approach enabled us to assess the algorithm's performance across a range of conditions and noise levels. The uncertainty analysis was performed through three synthetic models, each incorporating varying noise levels, regional effects, neighboring structure influences, and variable window lengths. The relative error and standard deviation for each model parameter were calculated and presented in Tables 2, 3, 4, 5, 6, 7, 8, 9 and 10; Figs. 7, 10, and 15.
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Third model, two neighboring structures, in case of a 10% noise level condition. Panel (a) exhibits the observed gravity anomaly profile (illustrated by black circles) alongside the computed optimal gravity response (depicted by red circles) employing the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using varying ‘s’ window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

<table><caption>**Table 9.** Results of parameter estimation for the two neighboring synthetic models under a noise-free condition.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s=1 km</th><th>s=2 km</th><th>s=3 km</th><th>s=4 km</th><th>s=5 km</th><th>s=6 km</th><th>s=7 km</th><th>s=8 km</th><th>s=9 km</th><th>s=10 km</th></tr></thead><tbody><tr><td>A<sub>1</sub> (mGal.km)</td><td>230</td><td>100-600</td><td>255</td><td>240</td><td>215</td><td>230</td><td>235</td><td>240</td><td>235</td><td>245</td><td>235</td><td>220</td><td>235 ± 11.547</td><td rowspan="8">2.15*10<sup>-6</sup></td></tr><tr><td>z<sub>1</sub> (km)</td><td>5</td><td>1-10</td><td>5.5</td><td>5.4</td><td>5.4</td><td>4.9</td><td>4.9</td><td>4.8</td><td>5.3</td><td>5.4</td><td>4.7</td><td>4.7</td><td>5.1 ± 0.327</td></tr><tr><td>x<sub>o1</sub> (km)</td><td>-20</td><td>-30-30</td><td>-9</td><td>-8</td><td>-11</td><td>-10</td><td>-10</td><td>-11</td><td>-8</td><td>-9</td><td>-11</td><td>-13</td><td>-10 ± 1.564</td></tr><tr><td>q<sub>1</sub></td><td>0.5</td><td>0-2</td><td>0.52</td><td>0.52</td><td>0.48</td><td>0.47</td><td>0.49</td><td>0.48</td><td>0.49</td><td>0.5</td><td>0.53</td><td>0.52</td><td>0.5 ± 0.021</td></tr><tr><td>m<sub>1</sub></td><td>0</td><td>0-2</td><td>0.02</td><td>0.01</td><td>0</td><td>0</td><td>0.01</td><td>0.02</td><td>0.02</td><td>0</td><td>0.02</td><td>0.01</td><td>0.011 ± 0.009</td></tr><tr><td>A<sub>2</sub> (mGal.km)</td><td>200</td><td>100-600</td><td>195</td><td>210</td><td>205</td><td>195</td><td>195</td><td>200</td><td>195</td><td>180</td><td>185</td><td>190</td><td>195 ± 8.819</td></tr><tr><td>z<sub>2</sub> (km)</td><td>3</td><td>1-10</td><td>2.8</td><td>2.7</td><td>3.1</td><td>3.2</td><td>3.1</td><td>3</td><td>2.8</td><td>2.7</td><td>2.8</td><td>2.8</td><td>2.9 ± 0.183</td></tr><tr><td>x<sub>o2</sub> (km)</td><td>25</td><td>-30-30</td><td>12</td><td>8</td><td>11</td><td>10</td><td>9</td><td>10</td><td>10</td><td>9</td><td>10</td><td>11</td><td>10 ± 1.155</td></tr><tr><td>q<sub>2</sub></td><td>0.5</td><td>0-2</td><td>0.45</td><td>0.48</td><td>0.49</td><td>0.5</td><td>0.52</td><td>0.54</td><td>0.5</td><td>0.52</td><td>0.52</td><td>0.48</td><td>0.5 ± 0.026</td></tr><tr><td>m<sub>2</sub></td><td>0</td><td>0-2</td><td>0.02</td><td>0.02</td><td>0</td><td>0.01</td><td>0</td><td>0.01</td><td>0.02</td><td>0.01</td><td>0</td><td>0.01</td><td>0.01 ± 0.008</td></tr></tbody></table>

In the first synthetic model, the algorithm was tested under three noise conditions: noise-free, 5% random noise, and 10% random noise. This model allowed us to assess the impact of different noise levels on parameter estimation accuracy. Figure 7 demonstrates the algorithm's robustness and accuracy across these noise levels, illustrating its ability to maintain reliable performance even in the presence of significant noise.

In the second synthetic model, the effect of a first-order regional anomaly was examined under both noise-free and 10% random noise conditions. Figure 10 showcases the algorithm's ability to accurately estimate parameters despite regional anomalies and noise, highlighting its effectiveness in isolating local anomalies from broader regional trends.

The third synthetic model focused on the influence of neighboring structures on parameter estimation. Two scenarios were tested: two separated structures and two neighboring structures, with both noise-free and
<table><caption>**Table 10.** Results of parameter estimation for the two neighboring synthetic model under a scenario with a 10% noise level.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">True model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s=1 km</th><th>s=2 km</th><th>s=3 km</th><th>s=4 km</th><th>s=5 km</th><th>s=6 km</th><th>s=7 km</th><th>s=8 km</th><th>s=9 km</th><th>s=10 km</th></tr></thead><tbody><tr><td>A<sub>1</sub> (mGal.km)</td><td>230</td><td>100-600</td><td>240</td><td>235</td><td>220</td><td>235</td><td>225</td><td>235</td><td>220</td><td>210</td><td>225</td><td>205</td><td>225 ± 11.547</td><td rowspan="8">1.16*10<sup>-5</sup></td></tr><tr><td>z<sub>1</sub> (km)</td><td>5</td><td>1-10</td><td>5.5</td><td>5.4</td><td>5.3</td><td>5.3</td><td>5.2</td><td>5.2</td><td>5.3</td><td>5.2</td><td>5.3</td><td>5.3</td><td>5.3 ± 0.094</td></tr><tr><td>x<sub>o1</sub> (km)</td><td>-20</td><td>-30-30</td><td>-13</td><td>-12</td><td>-12</td><td>-10</td><td>-10</td><td>-9</td><td>-11</td><td>-10</td><td>-11</td><td>-12</td><td>-11 ± 1.247</td></tr><tr><td>q<sub>1</sub></td><td>0.5</td><td>0-2</td><td>0.5</td><td>0.52</td><td>0.52</td><td>0.5</td><td>0.48</td><td>0.48</td><td>0.5</td><td>0.46</td><td>0.48</td><td>0.46</td><td>0.49 ± 0.022</td></tr><tr><td>m<sub>1</sub></td><td>0</td><td>0-2</td><td>0.01</td><td>0</td><td>0.02</td><td>0.02</td><td>0</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0</td><td>0.01</td><td>0.01 ± 0.008</td></tr><tr><td>A<sub>2</sub> (mGal.km)</td><td>200</td><td>100-600</td><td>220</td><td>215</td><td>200</td><td>210</td><td>210</td><td>215</td><td>200</td><td>195</td><td>215</td><td>220</td><td>210 ± 8.819</td></tr><tr><td>z<sub>2</sub> (km)</td><td>3</td><td>1-10</td><td>2.7</td><td>3.1</td><td>3.1</td><td>2.8</td><td>2.7</td><td>2.8</td><td>2.9</td><td>2.7</td><td>2.6</td><td>2.6</td><td>2.8 ± 0.183</td></tr><tr><td>x<sub>o2</sub> (km)</td><td>25</td><td>-30-30</td><td>11</td><td>8</td><td>12</td><td>8</td><td>9</td><td>11</td><td>10</td><td>12</td><td>11</td><td>8</td><td>10 ± 1.633</td></tr><tr><td>q<sub>2</sub></td><td>0.5</td><td>0-2</td><td>0.52</td><td>0.52</td><td>0.48</td><td>0.5</td><td>0.48</td><td>0.48</td><td>0.5</td><td>0.54</td><td>0.52</td><td>0.46</td><td>0.5 ± 0.025</td></tr><tr><td>m<sub>2</sub></td><td>0</td><td>0-2</td><td>0.01</td><td>0</td><td>0</td><td>0.01</td><td>0</td><td>0.01</td><td>0</td><td>0</td><td>0.01</td><td>0.01</td><td>0.005 ± 0.005</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** The robustness and accuracy of the parameter estimations for the third model, in case separated and neighboring structures without and with 10% random noise.</figcaption></figure>

10% random noise conditions. Figure 15 confirms the algorithm's stability and accuracy in handling complex geological settings with closely spaced structures. This model further validates the algorithm's robustness in diverse subsurface environments. The percentage error in parameter estimation was calculated, and the results are visually represented in the corresponding figures, reinforcing the algorithm's reliability and precision.

Overall, the close agreement between the recovered parameters and the original gravity response, even in the presence of noise and complex geological conditions, demonstrates the MBA's reliability and effectiveness. These results suggest that the MBA approach is well-suited for interpreting real-world gravity data impacted by various uncertainties.

## Real datasets
Mobrun Anomaly, Québec, Canada

The Mobrun gravity anomaly located approximately 34 km northeast of Rouyn-Noranda in Québec, Canada, exposes a varied geological landscape characterized by the presence of the Renault Formation³⁸,³⁹. Bordered by Destor rocks running from north to south, the Renault Formation displays a mix of large fragmental rhyolitic and
andesitic units, contributing to the area's geological intricacy³⁸,³⁹ (Figs. 16 and 17). In the footwall of the Mobrun area, the Copper Hill rhyolite unit predominates, followed by a complex sequence of andesitic-rhyolitic rocks, felsic pyroclastic materials, and a layer of massive rhyolite forming the hanging wall³⁹,⁴⁰. Beneath the primary complex lies brecciated rhyolite, adding another layer of geological complexity to the region³⁹,⁴⁰. Notably, the Mobrun deposit, a key point of geological interest in the area, contains two extensive sulfide lens complexes, with the primary complex hosting five significant ore bodies³⁹.

A 230 m Gravity profile is taken across the gravity map (Fig. 18)³⁸ and shown in Fig. 19a. Firstly, the SHG method was employed using different s values (for s = 1 to 10 m) (Fig. 19b) to remove the regional effect. Then, the MBA scheme was applied to estimate pivotal parameters characterizing the Quebec anomaly. The global and average NRMSE solutions of the objective function are shown in Fig. 19c and d, respectively. The parameters resulting from various window values and their averages are meticulously detailed in Table 11 for comprehensive evaluation. The estimated parameters are summarized as follows: A = 80 ± 2.31 mGal.m, z = 48 ± 1.16 m, x₀ = 0 ± 1.15 m, q = 1 ± 0.08, m = 1 ± 0.07, these parameters collectively characterize a horizontal cylinder. Furthermore, the objective function value corresponding to these estimations is calculated to be 0.00016. To contextualize our findings, a comparative analysis with prior research has been undertaken, revealing the following insights. The summarized comparisons with existing studies are presented in Table 12, juxtaposing our estimated parameters with those from Mehanee⁴⁰, Biswas³³, Singh and Biswas⁴¹, Essa et al.⁴², and Elhussein and Diab⁴³. These systematic comparisons illuminate the efficacy of our algorithm in accurately estimating parameters, thus augmenting our understanding of the Quebec anomaly's geological characteristics and contributing to the broader geological discourse.

## Phenaimata gravity anomaly, Gujarat, India

The Phenaimata gravity anomaly, situated in the Narmada-Tapti tectonic zone (NTTZ) in Gujarat, India, poses a fascinating geological puzzle where igneous complexity intersects with tectonic dynamics. Situated north of the Narmada River, the Phenaimata igneous complex stands as a testament to the intricate geological formations shaped by the region's tectonic evolution⁴⁴. Comprising a mix of plutonic and volcanic rock types, including basalt, layered gabbro, diorite, nepheline-syenite, lamprophyres, and granophyres, the Phenaimata plug showcases a dynamic interplay of magmatic processes⁴⁵,⁴⁶ (Figs. 20 and 21). The particular significance lies in the prevalence of orthopyroxene gabbro within the Phenaimata igneous complex, representing a crucial stage

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** The geological map of the Noranda property in Canada, adapted from³⁸. The open squares symbolize the mines currently in operation, while the closed squares denote the mines that have ceased production.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Geological Map of the Mobrun Property shows the Rock Units and Structural Features (After³⁸).</figcaption></figure>

in the magmatic evolution and differentiation of the intrusive body⁴⁴. Petrological modeling suggests that the formation of gabbroic rocks within the complex is intricately linked to magma accumulation in the crust, along with assimilation and fractional crystallization processes⁴⁶.

Geophysical investigations, including Bouguer gravity anomaly mapping, have unveiled the distinctive gravitational signature of the Phenaimata igneous exposure⁴⁷,⁴⁸. An elliptical-shaped anomaly closure to the north of the complex aligns with the orientation of the Narmada-Tapti tectonic zone, underscoring the tectonic control exerted on the regional gravity field⁴⁷,⁴⁸.

To unravel the subsurface architecture and gravitational anomalies accompanying the Phenaimata complex, a meticulous north-south BB' 31 km profile was extracted perpendicular to the closure anomaly in the Bouguer gravity map⁴⁷,⁴⁸ (Fig. 22). Through the judicious application of separation techniques, the BB' residual gravity anomaly profile was derived (Fig. 23a), providing invaluable insights into the geological underpinnings shaping the Phenaimata gravity anomaly⁴⁷,⁴⁸.

Firstly, the SHG method was employed using different s values (for s=1 to 10 km) (Fig. 23b) to remove the regional effect. Then, the MBA process was applied to estimate pivotal parameters (Table 13). The global and average NRMSE solutions of the objective function are shown in Fig. 23c and d, respectively. The estimated parameters are summarized as follows: A = 1095 ± 5.27 mGal.km², z = 5.45 ± 0.08 km, x₀ = 1 ± 1.41 km, q = 1.5 ± 0.06, and m = 1 ± 0.06, these parameters collectively characterize a spherical body. These estimations represent a significant advancement in our understanding of the Phenaimata anomaly's geological characteristics.

For comparison with previous studies, a comprehensive comparison was conducted, incorporating findings from Essa and Diab¹⁹ (Table 14). This systematic analysis underscores the efficacy and reliability of our algorithm in accurately estimating parameters, contributing to a deeper understanding of the Phenaimata anomaly's geological nuances. With these results, significant progress has been made in unraveling the geological complexities of the Phenaimata anomaly and enhancing our understanding of Gujarat's broader geological landscape.
<figure><img src="image_20.png" /><figcaption>Fig. 18. The Mobrun gravity anomaly map (modified and redrawn after³⁸). The line marked as AB represents the gravity anomaly profile (see Fig. 19) that is open to interpretation.</figcaption></figure>

## Camaguey chromite anomaly, Cuba

The chromite deposits found in the Camaguey District of Cuba present a fascinating geological mystery. Here, serpentinized dunite and peridotite rocks intersect with feldspathic lithologies, creating a unique geological setting that is ripe for exploration⁴⁹⁻⁵¹ (Fig. 24).

In a groundbreaking endeavor, Davis et al.⁴⁹ led an exploration campaign under the auspices of the United States Geological Survey. Their goal was to meticulously collect gravity data to uncover the mysterious subsurface structure of the chromite-bearing ore deposits in the Camaguey District. This detailed data collection initiative aimed to map out the structural framework and density anomalies associated with the chromite mineralization, providing crucial insights into the geological origins and economic potential of the deposits.

The resultant gravity map (Fig. 25)⁴⁹, meticulously crafted over the Camaguey chromite deposits, serves as a cartographic testament to the intricate interplay of gravitational forces and geological heterogeneities shaping subterranean realm⁴⁹. Figuratively etched onto this geological canvas are the spatial delineations and density gradients delineating the chromite-bearing lithologies, signifying the gravitational im-print of the underlying mineralization.

To augment the understanding of the subsurface structure, a gravity profile perpendicular to the chromite deposit's strike was meticulously extracted from the residual gravity map⁵⁰. The profile spanning a length of 89 m (Fig. 26a), this gravity profile serves as a veritable Rosetta Stone, unlocking the gravitational secrets harbored within the subterranean realm.

By applying our suggested procedure, which started by applying Second Horizontal Derivative using different s values (for s=1 to s=10 m) (Fig. 26b) Then, the Bat Algorithm optimization was applied to estimate pivotal parameters characterizing the Camaguey anomaly. The global and average NRMSE solutions of the objective function are shown in Fig. 26c and d, respectively. The estimated parameters, meticulously obtained through the application of various window values and their averages, are presented in Table 15 for comprehensive assessment. For the Camaguey anomaly, the estimated parameters are summarized as follows: A = 22.99 ± 0.19 mGal.m, z = 14.72 ± 0.22 m, x₀ = -1 ± 1.16 m, q = 1 ± 0.13, and m = 1 ± 0.06, these parameters collectively characterize a horizontal cylinder. These estimations offer valuable insights into the geological structure and characteristics of the Camaguey anomaly compared to the data from the drilled wells (Fig. 27). To facilitate a thorough comparison with existing studies, we have juxta-posed our estimated parameters with those reported in previous works. The summarized comparisons, inclusive of findings from Mehanee⁴⁰, Biswas³³, Ekinci et al.⁵⁰, Essa et al.⁴², and Elhussein and Diab⁴³, are meticulously presented in Table 16.
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** First field dataset, Mobrun gravity anomaly, interpretation. Panel (a) illustrates the observed gravity anomaly profile (marked by black circles), the computed optimal gravity response (indicated by blue circles) using the MBA algorithm technique, and the calculated optimal gravity response (depicted by purple asteroids) by$^{43}$. Panel (b) demonstrates the SHG anomalies using various 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

<table><caption>**Table 11.** Results of parameter estimation for the Mobrun anomaly, Quebec, Canada.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 m</th><th>s = 2 m</th><th>s = 3 m</th><th>s = 4 m</th><th>s = 5 m</th><th>s = 6 m</th><th>s = 7 m</th><th>s = 8 m</th><th>s = 9 m</th><th>s = 10 m</th></tr></thead><tbody><tr><td>A (mGal.m)</td><td>0–300</td><td>78</td><td>80</td><td>82</td><td>84</td><td>76</td><td>80</td><td>80</td><td>82</td><td>80</td><td>78</td><td>80 ± 2.31</td><td rowspan="5">0.00016</td></tr><tr><td>z (m)</td><td>10–100</td><td>47</td><td>48</td><td>46</td><td>47</td><td>49</td><td>48</td><td>50</td><td>48</td><td>48</td><td>49</td><td>48 ± 1.16</td></tr><tr><td>x<sub>o</sub> (m)</td><td>-10–10</td><td>-1</td><td>-2</td><td>0</td><td>0</td><td>1</td><td>1</td><td>2</td><td>0</td><td>0</td><td>-1</td><td>0 ± 1.15</td></tr><tr><td>q</td><td>0–2</td><td>1.05</td><td>1</td><td>0.9</td><td>0.95</td><td>1</td><td>1</td><td>1.15</td><td>0.9</td><td>0.95</td><td>1.1</td><td>1 ± 0.08</td></tr><tr><td>m</td><td>0–2</td><td>1</td><td>1</td><td>0.9</td><td>0.95</td><td>1.05</td><td>1.1</td><td>0.9</td><td>1</td><td>1.1</td><td>1</td><td>1 ± 0.07</td></tr></tbody></table>

<table><caption>**Table 12.** Comparative analysis of the retrieved model parameters for the Mobrun anomaly in Quebec, Canada, derived from our algorithm and those reported in previous studies.</caption><thead><tr><th>Model parameters</th><th>Mehanee <sup>40</sup></th><th>Biswas <sup>33</sup></th><th>Singh and Biswas <sup>41</sup></th><th>Essa et al. <sup>42</sup></th><th>Elhussein and Diab <sup>43</sup></th><th>Present Study</th></tr></thead><tbody><tr><td>A (mGal.m)</td><td>80</td><td>79.5</td><td></td><td>38.47</td><td>75.33</td><td>80 ± 2.31</td></tr><tr><td>z (m)</td><td>47</td><td>47.7</td><td>46.69</td><td>47</td><td>46</td><td>48 ± 1.16</td></tr><tr><td>xo (m)</td><td></td><td>2.5</td><td>2.37</td><td>2</td><td>2</td><td>0 ± 1.15</td></tr><tr><td>q</td><td>1</td><td>1</td><td>1</td><td>0.91</td><td>1</td><td>1 ± 0.08</td></tr><tr><td>m</td><td></td><td></td><td></td><td></td><td></td><td>1 ± 0.07</td></tr></tbody></table>

These systematic comparisons shed light on the efficacy and reliability of our algorithm in accurately estimating parameters, thus contributing to a deeper understanding of the Camaguey anomaly's geological intricacies.

Finally, our algorithm demonstrates excellent performance across a variety of synthetic scenarios with differing complexities, including variations in noise levels, regional influences, and multi-model effects. Its
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** The geological and tectonic features of the Deccan traps surrounding the Phenaimata igneous intrusion (After$^{45}$).</figcaption></figure>

efficacy is further validated through real-world applications, such as the Mobrun sulfide deposit in Canada, the Phenaimata igneous complex intrusion in India, and the Camaguey chromite deposits in Cuba.

The algorithm effectively balances the exploration and exploitation phases, achieving rapid convergence and adaptability to complex optimization challenges. This is largely due to the integration of the MBA and SHG techniques, which enhance the delineation of subsurface anomalies. By accentuating anomaly boundaries and minimizing the impact of regional gravity fields, the algorithm facilitates more precise and reliable parameter estimation for subsurface structures.

However, despite these strengths, the method has certain limitations. Its computational demands can be significant, particularly when applied to large datasets or highly complex geological settings. Furthermore, the algorithm's performance is sensitive to parameter configurations, requiring meticulous tuning to optimize its outcomes.

## Conclusions

Our proposed algorithm has demonstrated exceptional effectiveness in accurately estimating parameters related to subsurface geological bodies, both in real-world and synthetic scenarios. Through rigorous analysis and optimization, we have achieved precise estimations for crucial parameters, including amplitude coefficient, depth, horizontal location, shape factor, and shape change coefficient. The validation of our algorithm across diverse geological settings such as anomalies in Quebec (Canada), Phenaimata (India), and Camaguey (Cuba), as well as through synthetic models, consistently yielded accurate results, underscoring the robustness and reliability of our approach in handling complex geological data. We strongly recommend that researchers adopt our algorithm for geological investigations, as it provides a powerful tool for accurately estimating subsurface
<figure><img src="image_23.png" /><figcaption>**Fig. 21.** The geological features of the Deccan traps surrounding the Phenaimata igneous intrusion (After⁴⁵). The location of the gravity profile shown in the figure as BB'.</figcaption></figure>

geological parameters. Our approach outperforms existing methodologies in both accuracy and efficiency, incorporating advanced techniques like the second horizontal gradient (SHG) method and the metaheuristic Bat algorithm (MBA). Furthermore, the adaptability of our algorithm to various geological scenarios, noise levels, and data complexities emphasizes its versatility and applicability across a wide range of research domains. We encourage researchers to incorporate our algorithm into their workflows to enhance the accuracy and reliability of their geological analyses and interpretations. The adaptability and versatility of our algorithm highlight the potential for innovative technology to revolutionize research in deciphering subsurface geological structures with unprecedented precision and efficiency.
<figure><img src="image_24.png" /><figcaption>**Fig. 22.** This figure shows the original data of the profile BB' and a geological model beneath the gravity profile (After$^{45,52}$).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 23.** Second field dataset, Phenaimata gravity anomaly, interpretation. Panel (a) illustrates the observed gravity anomaly profile (marked by black circles) and the computed optimal gravity response (indicated by blue circles) using the MBA algorithm technique. Panel (b) demonstrates the SHG anomalies using various 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obi}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>
<table><caption>**Table 13.** Results of parameter estimation for the Phenaimata anomaly, Gujarat, India.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 km</th><th>s = 2 km</th><th>s = 3 km</th><th>s = 4 km</th><th>s = 5 km</th><th>s = 6 km</th><th>s = 7 km</th><th>s = 8 km</th><th>s = 9 km</th><th>s = 10 km</th></tr></thead><tbody><tr><td>A (mGal.km<sup>2</sup>)</td><td>800–1300</td><td>1090</td><td>1090</td><td>1105</td><td>1095</td><td>1090</td><td>1100</td><td>1090</td><td>1095</td><td>1095</td><td>1100</td><td>1095 ± 5.27</td><td rowspan="4">0.0000024</td></tr><tr><td>z (km)</td><td>1–10</td><td>5.4</td><td>5.45</td><td>5.5</td><td>5.5</td><td>5.6</td><td>5.3</td><td>5.45</td><td>5.5</td><td>5.4</td><td>5.4</td><td>5.45 ± 0.08</td></tr><tr><td>x<sub>o</sub> (km)</td><td>-10–10</td><td>2</td><td>1</td><td>2</td><td>-1</td><td>3</td><td>0</td><td>2</td><td>-1</td><td>0</td><td>2</td><td>1 ± 1.41</td></tr><tr><td>q</td><td>0–2</td><td>1.5</td><td>1.5</td><td>1.45</td><td>1.4</td><td>1.6</td><td>1.55</td><td>1.5</td><td>1.45</td><td>1.5</td><td>1.55</td><td>1.5 ± 0.06</td></tr><tr><td>m</td><td>0–2</td><td>0.95</td><td>0.95</td><td>0.95</td><td>1.05</td><td>1</td><td>0.95</td><td>1.05</td><td>1.05</td><td>0.95</td><td>1.1</td><td>1 ± 0.06</td><td></td></tr></tbody></table>

<table><caption>**Table 14.** Comparative analysis of the retrieved model parameters for the Phenaimata anomaly in Gujarat, India, obtained through our algorithm and those reported in previous studies.</caption><thead><tr><th>Model parameters</th><th>Essa and Diab <sup>19</sup></th><th>Present Study</th></tr></thead><tbody><tr><td>A (mGal.km<sup>2</sup>)</td><td>960</td><td>1095 ± 5.27</td></tr><tr><td>z (km)</td><td>5.3</td><td>5.45 ± 0.08</td></tr><tr><td>x<sub>o</sub> (km)</td><td>1</td><td>1 ± 1.41</td></tr><tr><td>q</td><td>1.5</td><td>1.5 ± 0.06</td></tr><tr><td>m</td><td>1</td><td>1 ± 0.06</td></tr></tbody></table>
<figure><img src="image_26.png" /><figcaption>**Fig. 24.** The geographical position and geological characteristics of the Camaguey District after⁵¹.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 25. The Camaguey gravity anomaly map redrawn and modified after⁴⁹. The line marked as AB represents the gravity anomaly profile (see Fig. 26) that is open to interpretation.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. 26.** Third field dataset, Camaguey gravity anomaly, interpretation. Panel (a) illustrates the observed gravity anomaly profile (marked by black circles), and the computed optimal gravity response (indicated by blue circles) using the MBA algorithm technique, and the calculated optimal gravity response (depicted by purple asteroids) by$^{43}$. Panel (b) demonstrates the SHG anomalies using various 's' window values. Panel (c) depicts the NRMSE of the global optimum solution ($F_{Obj}$) of the bats relative to the number of iterations. Panel (d) shows the average NRMSE of all the bats.</figcaption></figure>

<table><caption>**Table 15.** Results of parameter estimation for r the Camaguey anomaly, Cuba.</caption><thead><tr><th rowspan="2">Model parameters</th><th rowspan="2">Search range</th><th colspan="10">Estimated parameters</th><th rowspan="2">Avg. value ± Uncertainty</th><th rowspan="2">F<sub>obj</sub></th></tr><tr><th>s = 1 m</th><th>s = 2 m</th><th>s = 3 m</th><th>s = 4 m</th><th>s = 5 m</th><th>s = 6 m</th><th>s = 7 m</th><th>s = 8 m</th><th>s = 9 m</th><th>s = 10 m</th></tr></thead><tbody><tr><td>A (mGal.m)</td><td>0–50</td><td>3</td><td>2.7</td><td>3</td><td>2.8</td><td>3.1</td><td>3.1</td><td>3.1</td><td>3.2</td><td>3.2</td><td>2.7</td><td>2.99 ± 0.19</td><td rowspan="5">0.00012</td></tr><tr><td>z (m)</td><td>10–20</td><td>14.5</td><td>14.5</td><td>14.7</td><td>14.8</td><td>14.7</td><td>15</td><td>14.5</td><td>15</td><td>15</td><td>14.5</td><td>14.72 ± 0.22</td></tr><tr><td>x<sub>o</sub> (m)</td><td>-10–-10</td><td>0</td><td>-1</td><td>-1</td><td>-1</td><td>-2</td><td>1</td><td>0</td><td>-2</td><td>-3</td><td>-1</td><td>-1 ± 1.16</td></tr><tr><td>q</td><td>0–2</td><td>0.9</td><td>1.1</td><td>1.1</td><td>1.2</td><td>0.8</td><td>0.9</td><td>0.9</td><td>1</td><td>1</td><td>1.1</td><td>1 ± 0.13</td></tr><tr><td>m</td><td>0–2</td><td>0.95</td><td>1</td><td>0.95</td><td>1</td><td>1</td><td>1.05</td><td>1.1</td><td>0.9</td><td>1</td><td>1.05</td><td>1 ± 0.06</td></tr></tbody></table>
<figure><img src="image_29.png" /><figcaption>Fig. 27. shows the drilling data of the Camaguey chromite ore body, Cuba. (redrawn from Fig. 4;⁴⁹).</figcaption></figure>

<table><caption>**Table 16.** Comparative analysis of the retrieved model parameters for the Camaguey anomaly, Cuba, derived from our algorithm and those reported in previous studies.</caption><thead><tr><td>Model parameters</td><td>Mehanee <sup>40</sup></td><td>Biswas <sup>33</sup></td><td>Ekinci et al. <sup>50</sup></td><td>Essa et al. <sup>42</sup></td><td>Elhussein and Diab <sup>43</sup></td><td>Present Study</td></tr></thead><tbody><tr><td>A (mGal.m)</td><td>3</td><td>3.5</td><td>111.34 (mGal.m<sup>2</sup>)</td><td>1.86</td><td>3.03</td><td>2.99 ± 0.19</td></tr><tr><td>z (m)</td><td>16</td><td>16.2</td><td>23.23</td><td>16</td><td>15</td><td>14.72 ± 0.22</td></tr><tr><td>xo (m)</td><td></td><td>-1.8</td><td>58.73</td><td>0</td><td>-2</td><td>-1 ± 1.16</td></tr><tr><td>q</td><td>1</td><td>1</td><td>1.5</td><td>0.9</td><td>1</td><td>1 ± 0.13</td></tr><tr><td>m</td><td></td><td></td><td></td><td></td><td></td><td>1 ± 0.06</td></tr></tbody></table>

## Data availability

The authors declare that the data is available upon request from the corresponding author.

Received: 16 August 2024; Accepted: 28 January 2025

Published online: 10 February 2025

## References

1. Anderson, N. L., Essa, K. S. & Elhussein, M. A comparison study using particle swarm optimization inversion algorithm for gravity anomaly interpretation due to a 2D vertical fault structure. *J. Appl. Geophys.* **179**, 104120 (2020).

2. Ombati, D., Githiri, J., K'Orowe, M. & Nyakundi, E. Delineation of Subsurface Structures Using Gravity Data of the Shallow Offshore, Lamu Basin, Kenya. *Int. J. Geophys.* **3024977** (2022).

3. Essa, K. S. & Diab, Z. E. Exploring Fault Plane geometry through metaheuristic Bat Algorithm (MBA) analysis of potential Field Data: Environmental and Engineering Applications. *Rock. Mech. Rock. Eng.* https://doi.org/10.1007/s00603-024-04198-6 (2024).
4. Elhussein, M., Abo-Ezz, E. R., Gomaa, O. A., Géraud, Y. & Essa, K. S. Two-sided fault parameter assessment through the utilization of a particle swarm optimizer to different gravity horizontal gradients-orders with application to various cases studies. *Environ. Earth Sci.* **82**, 502 (2023).

5. Sadraeifar, B. & Abedi, M. A comparative study of sparse and Tikhonov regularization methods in gravity inversion: A case study of manganese deposit in Iran. *Int. J. Min. Geo-Eng.* **58**, 161–169 (2024).

6. Warnana, D. D., Sungkono, S. & Essa, K. S. Efficient and robust estimation of various Ore and Mineral Model parameters from residual gravity anomalies using the dual classification learning Rao Algorithm. *Arab. J. Sci. Eng.* https://doi.org/10.1007/s13369-024-09774-0 (2024).

7. Abou Zaki, N., Haghighi, T., Rossi, A. M., Tourian, P. J., Kløve, B. & M. & Monitoring Groundwater Storage Depletion using gravity recovery and climate experiment (GRACE) data in Bakhtegan Catchment. *Iran. Water.* **11**, 1456 (2019).

8. Essa, K. S. & Géraud, Y. Parameters estimation from the gravity anomaly caused by the two-dimensional horizontal thin sheet applying the global particle swarm algorithm. *J. Petrol. Sci. Eng.* **193**, 107421 (2020).

9. Eshagh, M. The Earth's Gravity Field Role in Geodesy and Large-Scale Geophysics. Geodetic Sciences Theory, Applications and Recent Developments. IntechOpen. (2021).

10. Essa, K. S., Mehanee, S. A. & Elhussein, M. Gravity data interpretation by a two-sided fault-like geologic structure using the global particle swarm technique. *Phys. Earth Planet. Inter.* **311**, 106631 (2021).

11. Guglielmetti, L. & Moscariello, A. On the use of gravity data in delineating geologic features of interest for geothermal exploration in the Geneva Basin (Switzerland): prospects and limitations. *Swiss J. Geosci.* **114**, 15 (2021).

12. Mohamed, A., Al Deep, M., Abdelrahman, K. & Abdelrady, A. Geometry of the Magma Chamber and Curie Point depth beneath Hawaii Island: inferences from magnetic and gravity data. *Front. Earth Sci.* **10**, 847984 (2022).

13. Velicogna, I. & Wahr, J. Time-variable gravity observations of ice sheet mass balance: Precision and limitations of the GRACE satellite data. *Geophys. Res. Lett.* **40**, 3055–3063 (2023).

14. Hirt, C. Gravity forward modelling. In (ed Graffarend, E.) Encyclopaedia of Geodesy. Curtin University, Perth, WA, Australia (2015).

15. Zhou, S., Wei, Y., Lu, P., Jiao, J. & Jia, H. Deep-learning gravity inversion method with depth-weighting constraints and its application in Geothermal Exploration. *Remote Sens.* **16** (23), 4467. https://doi.org/10.3390/rs16234467 (2024).

16. Essa, K. S. Gravity data interpretation using the s-curves method. *J. Geophys. Eng.* **4** (2), 204–213. https://doi.org/10.1088/1742-2132/4/2/2009 (2007).

17. Witter, J. B., Siler, D. L., Faulds, J. E. & Hinz, N. H. 3D geophysical inversion modeling of gravity data to test the 3D geologic model of the bradys geothermal area, Nevada, USA. *Geotherm. Energy* **4**(14) (2016).

18. Liu, Y. et al. Gravity predictions in data-missing areas using machine learning methods. *Remote Sens.* **16** (22), 4173 (2024).

19. Essa, K. S. & Diab, Z. E. Source parameters estimation from gravity data using Bat algorithm with application to geothermal and volcanic activity studies. *Int. J. Environ. Sci. Technol.* **20**, 4167–4187. https://doi.org/10.1007/s13762-022-04263-z (2023).

20. Devika, G. & Karegowda, A. G. Bio-inspired optimization: Algorithm, analysis and scope of application. In Swarm Intelligence - Recent Advances and Current Applications (2023).

21. Dubey, C. P., Pandey, L. & Rajalakshmi, K. V. Gradient descent fusion for gravity and magnetic data. *J. Earth Syst. Sci.* **133**, 124 (2024).

22. Ifiemi, E., Sodiki, J. & Nkoi, B. Analysis of a water distribution network by Newton-Raphson multivariable method: A case of negligible minor losses. *Innovative Syst. Des. Eng.* **11** (2), 15–28 (2020).

23. Biswas, A. Interpretation of gravity and magnetic anomaly over thin sheet-type structure using very fast simulated annealing global optimization technique. *Model. Earth Syst. Environ.* **2**, 30 (2016).

24. Niu, T., Zhang, G., Zhang, M. & Zhang, G. Joint inversion of gravity and gravity gradient data using smoothed L0 norm regularization algorithm with sensitivity matrix compression. *Front. Earth Sci.* **11**, 1283238 (2023).

25. Su, K. et al. Hunger games search for the elucidation of gravity anomalies with application to geothermal energy investigations and volcanic activity studies. *Open. Geosci.* **16** (1), 20220641 (2024).

26. Ai, H. et al. Y. Hunger games search optimization for the inversion of gravity anomalies of active mud diapir from SW Taiwan using inclined anticlinal source approximation. *J. Appl. Geophys.* **227**, 105443 (2024).

27. Vashist, A. & Mukhopadhyay, S. Stochastic inversion of gravity and magnetic data to build subsurface models using genetic algorithm, differential evolution, grey wolf optimizer, and particle swarm optimizer. *Stanf. Geotherm. Workshop* (2023).

28. Essa, K. S. & Munschy, M. Gravity data interpretation using the particle swarm optimisation method with application to mineral exploration. *J. Earth Syst. Sci.* **128**, 123 (2019).

29. Akila, M. & Christe, S. Multiobjective binary Bat Algorithm with greedy crossover for selecting the most promising features in classification problems. *Neural Comput. Appl.* **34** (12), 9876–9890 (2022).

30. Shehab, M., Abuhaija, B. et al. A comprehensive review of Bat inspired Algorithm: variants, applications, and hybridization. *Arch. Comput. Methods Eng.* **30**, 765–797 (2023).

31. Asfahani, J. & Tlas, M. An automatic method of direct interpretation of residual gravity anomaly profiles due to spheres and cylinders. *Pure appl. Geophys.* **165**, 981–994 (2008).

32. Essa, K. S. New fast least-squares algorithm for estimating the best-fitting parameters of some geometric-structures from gravity anomalies. *J. Adv. Res.* **5**, 57–65 (2014).

33. Biswas, A. Interpretation of residual gravity anomaly caused by simple shaped bodies using very fast simulated annealing global optimization. *Geosci. Front.* **6**, 875–893 (2015).

34. Telford, W. M., Geldart, L. P. & Sheriff, R. E. *Applied Geophysics* (Cambridge University Press, 1990).

35. Asfahani, J. & Tlas, M. The effects of the different window functions on the analytic signal and the horizontal gradient magnitude of potential field anomalies: A case study from Northern Lebanon. *J. Appl. Geophys.* **113**, 88–95 (2015).

36. Essa, K. S. & Elhussein, M. A new approach for the interpretation of magnetic data by a 2-D dipping dike. *J. Appl. Geophys.* **136**, 431–443 (2017).

37. Yang, X. S. A new metaheuristic bat-inspired algorithm. In: (eds Gonzalez, J., Pelta, D., Cruz, C., Terrazas, G. & Krasnogor, N.) Nature Inspired Cooperative Strategies. Springer, Berlin. (2010).

38. Barrett, T. J., Cattalani, S., Hoy, L., Riopel, J. & Lafleur, P. J. Massive sulfide deposits of the Noranda area, Quebec. IV. The Mobrun mine. *Can. J. Earth Sci.* **29**, 1349–1374 (1992).

39. Caumartin, C. et al. Volcanic stratigraphy and structure of the Mobrun mine. In the northwestern Quebec polymetallic belt. Edited by M. Rive, P.; Verpaelst, Y.; Gagnon, Canadian Institute of Mining and Metallurgy, 43, 133–142 (1990).

40. Mehanee, S. A. Accurate and efficient regularised inversion approach for the interpretation of isolated gravity anomalies. *Pure Appl. Geophys.* **171** (8), 1897–1937 (2014).

41. Singh, A. & Biswas, A. Application of global particle swarm optimization for inversion of residual gravity anomalies over geological bodies with idealized geometries. *Nat. Resour. Res.* **25**, 297–314 (2016).

42. Essa, K. S., Mehanee, S. A., Soliman, K. S. & Diab, Z. E. Gravity profile interpretation using the R-parameter imaging technique with application to ore exploration. *Ore Geol. Rev.* **126**, 103695 (2020).

43. Elhussein, M. & Diab, Z. E. Gravity Data Imaging Using Local Wavenumber-Based Algorithm: Sustainable Development Cases Studies. *Nat. Resour. Res.* **32**, 171–193 (2023).
44. Hari, K. R., Chalapathi Rao, N. V. & Swarnkar, V. Petrogenesis of gabbro and orthopyroxene gabbro from the Phenai mata igneous complex, Deccan volcanic province: products of concurrent assimilation and fractional crystallization. *J. Geol. Soc. India.* **78**, 501–509 (2011).

45. Sukeshwala, R. N. & Sethna, S. F. Layered gabbro of composite plug of Phenaimata, Gujarat State. *J. Geol. Soc. India.* **10**, 177–187 (1969).

46. Hari, K. R., Chalapathi Rao, N. V., Swarnkar, V. & Hou, G. Alkali feldspar syenites with shoshonitic afnities from Chhotaudepur area: implication for mantle metasomatism in the Deccan large igneous province. *Geosci. Front.* **5**, 261–276 (2014).

47. Basu, A. R., Renne, P. R., Dasgupta, D. K., Teichmann, F. & Poreda, R. J. Early & late alkali igneous pulses and a high-3He plume origin for the deccan food basalts. *Science* **261**, 902–906 (1993).

48. Poornachandra Rao, G. V. S., Mallikharjun Rao, J. & Jaya Prasanna Lakshmi, K. Palaeomagnetic study of the alkaline rocks associated with the Deccan traps of north western India. *Indian J. Geochem.* **19**, 19–32 (2004).

49. Davis, W. E., Jackson, W. H. & Richter, D. H. Gravity prospecting for chromite deposits in Camaguey Province. *Cuba Geophys.* **22**, 848–869 (1957).

50. Ekinci, Y. L., Balkaya, C., Gokturkler, G. & Turan, S. Model parameter estimations from residual gravity anomalies due to simple-shaped sources using Differential Evolution Algorithm. *J. Appl. Geophys.* **129**, 133–147 (2016).

51. Santana, M. M. U. et al. The La Unio'n au ± Custabled prospect, Camaguey District, Cuba: fluid inclusion and stable isotope evidence for ore-forming processes. *Min. Deposita.* **46** (1), 91–104 (2011).

52. Singh, B., Prabhakara Rao, M. R. K., Prajapati, S. K. & Swarnapriya, C. Combined gravity and magnetic modeling over Pavagadh and Phenaimata igneous complexes, Gujarat, India: inference on emplacement history of Deccan Volcanism. *J. Asian Earth Sci.* **80**, 119–133 (2014).

## Acknowledgements

The authors would like to thank Prof. Shreyas Deshpande, Editorial Board Member, and the three capable expert reviewers for their keen interest, valuable comments on the manuscript, and improvements to this work.

## Author contributions

All authors: Writing the original manuscript Draft, Conceptualization, Methodology, Figures and Tables preparations. All authors have read and agreed to the published version of the manuscript.

## Funding

This research received no external funding.

## Declarations

### Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to K.S.E.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025
Terms and Conditions

Springer Nature journal content, brought to you courtesy of Springer Nature Customer Service Center GmbH ("Springer Nature").

Springer Nature supports a reasonable amount of sharing of research papers by authors, subscribers and authorised users ("Users"), for small-scale personal, non-commercial use provided that all copyright, trade and service marks and other proprietary notices are maintained. By accessing, sharing, receiving or otherwise using the Springer Nature journal content you agree to these terms of use ("Terms"). For these purposes, Springer Nature considers academic use (by researchers and students) to be non-commercial.

These Terms are supplementary and will apply in addition to any applicable website terms and conditions, a relevant site licence or a personal subscription. These Terms will prevail over any conflict or ambiguity with regards to the relevant terms, a site licence or a personal subscription (to the extent of the conflict or ambiguity only). For Creative Commons-licensed articles, the terms of the Creative Commons license used will apply.

We collect and use personal data to provide access to the Springer Nature journal content. We may also use these personal data internally within ResearchGate and Springer Nature and as agreed share it, in an anonymised way, for purposes of tracking, analysis and reporting. We will not otherwise disclose your personal data outside the ResearchGate or the Springer Nature group of companies unless we have your permission as detailed in the Privacy Policy.

While Users may use the Springer Nature journal content for small scale, personal non-commercial use, it is important to note that Users may not:

1. use such content for the purpose of providing other users with access on a regular or large scale basis or as a means to circumvent access control;

2. use such content where to do so would be considered a criminal or statutory offence in any jurisdiction, or gives rise to civil liability, or is otherwise unlawful;

3. falsely or misleadingly imply or suggest endorsement, approval, sponsorship, or association unless explicitly agreed to by Springer Nature in writing;

4. use bots or other automated methods to access the content or redirect messages

5. override any security feature or exclusionary protocol; or

6. share the content in order to create substitute for Springer Nature products or services or a systematic database of Springer Nature journal content.

In line with the restriction against commercial use, Springer Nature does not permit the creation of a product or service that creates revenue, royalties, rent or income from our content or its inclusion as part of a paid for service or for other commercial gain. Springer Nature journal content cannot be used for inter-library loans and librarians may not upload Springer Nature journal content on a large scale into their, or any other, institutional repository.

These terms of use are reviewed regularly and may be amended at any time. Springer Nature is not obligated to publish any information or content on this website and may remove it or features or functionality at our sole discretion, at any time with or without notice. Springer Nature may revoke this licence to you at any time and remove access to any copies of the Springer Nature journal content which have been saved.

To the fullest extent permitted by law, Springer Nature makes no warranties, representations or guarantees to Users, either express or implied with respect to the Springer nature journal content and all parties disclaim and waive any implied warranties or warranties imposed by law, including merchantability or fitness for any particular purpose.

Please note that these rights do not automatically extend to content, data or other material published by Springer Nature that may be licensed from third parties.

If you would like to use or distribute our Springer Nature journal content to a wider audience or on a regular basis or in any other manner not expressly permitted by these Terms, please contact Springer Nature at

onlineservice@springernature.com