The geophysical inverse problem seeks to determine a model vector m ∈ ℝM from an observation vector d ∈ ℝN I L O ℝM ℝN denotes the forward modeling operator, which encodes the physics governing the relationship between subsurface

a Gaussian distribution with covariance matrix Cd. For seismic problems, G may involve the solution of the acoustic or elastic wave equation via finite-difference or finite-element methods; for gravity problems, G involves the evaluation of Newton's gravitational integral; and for Electromagnetic (EM) problems, G requires the solution of Maxwell's equations in heterogeneous media [29].

The inversion is typically cast as an optimization problem in which the objective function (also termed the cost function or fitness function) is minimized. The general form of the objective function is:

Φ(m) = Φd(m) + λ Φm(m), (1)

(m) = (d IT C 1 (d I I I

d obs d obs m

m )T C 1 I ) is the model regularization term that penalizes deviation from a reference model m ,

ref m ref ref

weighted by the model covariance matrix Cm. The trade- lance between fitting the data and maintaining model simplicity or geological plausibility [3], [4].

In the metaheuristic framework, this minimization is performed without computing gradients of respect to m. Instead, a population of candidate solutions (or a single solution in the case of trajectory-based methods such as SA) iteratively explores the model space through stochastic operators inspired by natural phenomena. The key advantage is that these methods require only the evaluation of the forward problem I I C-free" optimizers that can handle discontinuous, non-differentiable, and highly multimodal objective landscapes.

Metaheuristic algorithms applicable to geophysical problems can be broadly classified into four categories based on their inspiration source and search mechanism: evolutionary algorithms, swarm intelligence algorithms, physics-based algorithms, and hybrid/other nature-inspired algorithms. Table 1 provides a systematic classification of the principal algorithms that have been applied in geophysical contexts, along with their key characteristics relevant to inverse problem solving.

cooling schedule, handles high dimensionality; to cooling schedule; single-perturbation size robust for rugged landscapes solution trajectory GA Evolutionary 1975 Population size, Good global exploration; Premature convergence crossover rate, parallelizable; handles mixed risk; many control mutation rate, variable types parameters; selection method computationally expensive DE Evolutionary 1997 Scale factor F, Excellent for continuous Performance sensitive to F crossover rate CR, parameter spaces; few control and CR; may struggle with population size NP parameters; robust discrete parameters convergence PSO Swarm 1995 Inertia weight w, Fast convergence; simple Premature convergence in intelligence cognitive c₁, social implementation; excellent for multimodal spaces; poor c₂, population size continuous problems; few exploration in high parameters dimensions ACO Swarm 1992 Pheromone Naturally handles Not well-suited for intelligence evaporation rate, discrete/combinatorial continuous inversions; number of ants, problems; good for survey convergence speed issues heuristic weight design ABC Swarm 2005 Colony size, limit Balanced exploration- Slow convergence for high-intelligence parameter, number exploitation; minimal dimensional problems; of food sources parameters; robust for surface limited exploitation wave inversion capability CS Swarm 2009 Lévy flights enable efficient Limited exploitation; intelligence discovery exploration; simple; parameter sensitivity for probability pa, competitive performance step size population size BA Swarm 2010 Frequency range, Automatic balance between Convergence issues in high intelligence loudness A, pulse exploration and exploitation; dimensions; limited emission rate r echolocation-inspired local diversity maintenance search GWO Swarm 2014 Parameter vector a Parameter-free (no tuning May stagnate in complex intelligence (linearly beyond population size); good landscapes; limited diversity decreasing), convergence balance for large model spaces population size WOA Swarm 2016 Spiral shape Spiral updating enables fine- Slow convergence rate; may intelligence constant b, grained search; few parameters; miss global optimum in probability switch p good for 2D inversions high-dimensional spaces HHO Swarm 2019 Escaping energy E, Dynamic exploration- Relatively new; limited intelligence jump strength J, exploitation transition; geophysical validation; population size competitive with DE and PSO scalability unproven MFO Swarm 2015 Number of flames, Logarithmic spiral Premature convergence intelligence spiral constant b convergence; adaptive flame risk; limited exploration in reduction enhances later iterations exploitation SCA Math-based 2016 Constant a, random Simple trigonometric Limited performance on parameters r₁–r₄ operators; smooth transition complex multimodal from exploration to functions; fewer exploitation geophysical applications

The design of the objective function is critical to the success of any metaheuristic-based geophysical inversion, as it directly determines the topology of the search landscape and, consequently, the algorithm's ability to locate the global optimum. In geophysical practice, several norm-based formulations are employed. The L2- norm (least- d i (dobs,i calc,i)², is the most widely used due to its statistical optimality under Gaussian noise assumptions and its smooth, differentiable landscape. However, the L2-norm is highly sensitive to outliers and non-Gaussian noise, which are common in field geophysical data. The L1-norm d i |dobs,i calc,i|, provides greater robustness to outliers but introduces a non-differentiable objective landscape that is naturally suited to derivative-free metaheuristic optimization [30].

Hybrid norm formulations, such as the Huber norm, combine the sensitivity of L2 near the minimum with the robustness of L1 for large residuals. Additionally, geophysical objective functions may incorporate correlation-based measures (e.g., normalized cross-correlation for waveform matching in seismic inversion), envelope-based misfits for Full Waveform Inversion (FWI) that mitigate cycle-skipping issues, and multi-objective formulations that simultaneously minimize data misfit and geological complexity. For metaheuristic algorithms, the objective function evaluation requires a complete forward model computation for each candidate solution, making the computational cost of G(m) the primary bottleneck in the optimization process. This cost ranges from milliseconds for analytical gravity formulas to hours for three-dimensional elastic wave propagation simulations, fundamentally constraining the choice of algorithm and population size.

Regularization is essential in geophysical inversion to stabilize solutions and incorporate a priori geological information. In the metaheuristic framework, regularization can be implemented through several mechanisms. The most straightforward approach is the addition of a penalty term to the objective function, as described in Section 2.1, which penalizes model roughness (smoothness regularization), deviation from a reference model (Tikhonov regularization), or model complexity (minimum-structure regularization). The Occam approach, which seeks the smoothest model that fits the data to a prescribed tolerance, has been successfully adapted for metaheuristic inversions [4], [31].

A unique advantage of metaheuristic algorithms is their ability to naturally incorporate bound constraints by restricting the search space to physically plausible parameter ranges. Geological constraints such as layer ordering, minimum/maximum thickness, velocity-density relationships, and lithological associations can be enforced through constraint repair operators or penalty functions without modifying the underlying algorithm structure. This flexibility is particularly valuable in geophysics, where prior geological knowledge can dramatically reduce the non-uniqueness of the inverse problem. Additionally, metaheuristic populations provide an implicit form of uncertainty quantification: the ensemble of near-optimal solutions sampled during the optimization process approximates the posterior probability distribution of model parameters, enabling Bayesian-like uncertainty analysis without the formal machinery of Markov Chain Monte Carlo (MCMC) methods [32], [33].

Seismic methods represent the most widely used geophysical technique for subsurface exploration, and consequently, the largest body of metaheuristic applications in geophysics is found in seismic data processing and inversion. The review identified 64 publications (34.4% of the total) addressing seismic applications, spanning velocity model building, acoustic and Elastic Impedance (EI) inversion, Amplitude Variation with Offset (AVO) inversion, FWI, surface wave dispersion curve inversion, and seismic tomography.

FWI: FWI represents the most computationally demanding application of metaheuristics in seismology, as each forward model evaluation requires the numerical solution of the wave equation over the entire computational domain. Early work by Sen and Stoffa [23] demonstrated that VFSA could successfully recover one-dimensional velocity profiles from synthetic seismograms, establishing the feasibility of stochastic seismic inversion. Subsequent developments extended GA-based approaches to two-dimensional acoustic FWI, where populations of velocity models were evolved to minimize waveform misfits [35]. The PSO algorithm has been applied to FWI by multiple research groups, with notable success in mitigating the cycle-skipping problem that plagues gradient-based FWI when the initial model is far from the true solution [36]. Hybrid approaches combining DE with Gauss-Newton local optimization have shown particular promise, achieving convergence to accurate velocity models with 40–60% fewer forward simulations compared to standalone metaheuristic approaches.

Acoustic and EI inversion: post-stack and pre-stack seismic inversion for Acoustic Impedance (AI) and EI has been extensively addressed using metaheuristic algorithms. Kant et al. [37] presented a comprehensive comparison of GA and PSO for reservoir characterization through seismic inversion, demonstrating that PSO achieved a fitness error of 0.25 after 400 iterations compared to 0.88 for GA, with PSO also requiring 47% less computation time. The study applied both methods to synthetic data and real data from the Blackfoot field in Canada, successfully delineating a high-porosity reservoir zone (>20%) characterized by low AI (6000–8500 m/s·g/cc) in the 1040–1065 ms time range. DE has been applied to pre-stack simultaneous inversion for P-wave velocity, S-wave velocity, and density, exploiting AVO information to constrain elastic properties [38]. ABC algorithms have been employed for sparse-spike inversion with encouraging results in terms of resolution enhancement [39].

Surface wave dispersion inversion: the inversion of surface wave dispersion curves for near-surface shear-wave velocity profiling represents one of the most successful applications of metaheuristics in geophysics. The problem is inherently multimodal, with dispersion curves that can be fit equally well by substantially different velocity profiles, particularly when higher modes are present. PSO [40], GA [41], and SA [42] have all been applied successfully. More recently, GWO and WOA have been benchmarked against PSO for Rayleigh wave inversion, with GWO showing competitive accuracy and faster convergence for shallow profiles [43].

Seismic tomography: travel-time tomography and attenuation tomography have been addressed using metaheuristic approaches, primarily for cases where the ray geometry creates highly nonlinear relationships between velocity and travel time. Boschetti et al. [44] demonstrated GA-based seismic tomography, while Sambridge [45] introduced the Neighborhood Algorithm (NA) as a stochastic direct-search method for seismic tomographic inversion. Recent applications have employed DE and PSO for cross-borehole tomography with promising results for monitoring CO₂ sequestration sites [46].

Fig. 1 illustrates the convergence characteristics of four widely used metaheuristic algorithms, namely PSO, Genetic Algorithm (GA), DE, and VFSA, when applied to a benchmark two-dimensional AI inversion problem based on the Marmousi model. The horizontal axis represents the number of forward model evaluations ranging from 0 to 50,000, while the vertical axis shows the normalized logarithmic data misfit. PSO demonstrates the fastest initial convergence, achieving approximately 90% misfit reduction within the first 8,000 evaluations, indicating strong early-stage exploitation capability. In contrast, DE exhibits slower initial progress but achieves the lowest terminal misfit after approximately 30,000 evaluations, reflecting superior long-term global search performance and solution refinement capability. GA displays a characteristic step-wise convergence pattern caused by its generational selection and crossover operations, leading to intermittent improvements in solution quality. VFSA shows the smoothest and most monotonic convergence trajectory among all methods, although its convergence rate remains slower throughout the optimization process. Nevertheless, VFSA ultimately achieves highly competitive final solution quality, demonstrating its robustness for complex multimodal inverse problems. Overall, the figure highlights the trade-off between rapid convergence and final optimization accuracy across different metaheuristic paradigms in seismic inversion applications.

Potential field methods gravity and magnetics represent classical geophysical techniques that have benefited significantly from metaheuristic optimization. The inherent non-uniqueness of potential field inversion, arising from the fundamental ambiguity theorem (multiple source distributions can produce identical fields), makes these problems particularly well-suited to stochastic global optimization approaches. The review identified 42 publications (22.6%) addressing gravity and magnetic inversions.

Gravity anomaly inversion: the inversion of gravity anomalies for subsurface density structure has been approached using virtually every major metaheuristic algorithm. Shaw and Srivastava [26] provided one of the earliest comprehensive demonstrations of PSO for gravity inversion, estimating the geometry and density contrast of simple geological bodies (spheres, cylinders, and thin sheets) from residual gravity anomalies. Their work showed that PSO could reliably recover shape parameters to within 3–5% of true values from noise-free synthetic data and 8–12% from field data. The VFSA method was extensively used by Sen and Stoffa [47] for gravity inversion, demonstrating robust performance for estimating the depth, density contrast, and shape factor of buried geological structures.

Depth-to-basement estimation: the estimation of sedimentary basin depth from gravity data is a practically important problem that has been addressed by multiple metaheuristic approaches. DE-based inversions have been applied to determine basin geometry by discretizing the sediment-basement interface into a series of juxtaposed prisms, with DE optimizing the depths of individual prisms to minimize the misfit between observed and calculated Bouguer anomalies [48]. PSO variants with inertia weight adaptation have shown particular effectiveness for multi-body gravity inversion, where multiple geological sources contribute to the composite anomaly [49].

Magnetic susceptibility estimation: the inversion of magnetic anomalies for susceptibility distribution shares mathematical similarities with gravity inversion but introduces additional complexity through the vectorial nature of magnetization. Biswas and Acharya [50] demonstrated the application of ABC for magnetic anomaly interpretation, achieving depth estimation errors below 5% for synthetic models. The study of Vashisth et al. [51] at Stanford University presented a comprehensive comparison of DE, GA, PSO, and GWO for stochastic inversion of gravity and magnetic data to build subsurface geological fault models, finding that DE and PSO outperformed GA and GWO in terms of both convergence speed and solution accuracy for the multi-parameter fault geometry estimation problem.

EM geophysical methods, including Magnetotellurics (MT), Controlled-Source Electromagnetics (CSEM), Ground-Penetrating Radar (GPR), and Electrical Resistivity Tomography (ERT), provide crucial information about subsurface resistivity structure. The inverse problems associated with these methods are highly nonlinear and computationally intensive, as each forward model evaluation requires the numerical solution of Maxwell's equations in heterogeneous media. The review identified 31 publications (16.7%) focused on EM and resistivity inversions.

Magnetotelluric Inversion: the magnetotelluric method measures natural EM field variations to estimate Earth's resistivity structure from the surface to depths exceeding 100 km. The pioneering work of Dosso and Oldenburg [52] established the feasibility of SA-based MT inversion, applying SA to construct extremal conductivity models that bound the range of acceptable solutions. This work demonstrated that SA-based appraisal provided bounds comparable to those obtained via linearized inversion, with the critical advantage that SA results are not restricted to models near a reference solution. Subsequent developments applied VFSA to one-dimensional MT inversion [53], demonstrating that VFSA could efficiently sample the model space and recover resistivity profiles consistent with both synthetic and field data from Sundar Pahari, India. Comparative studies showed that PSO and GA provided results comparable to VFSA for 1D MT inversion but with faster convergence for shallow resistivity structures.

ERT: ERT inversion using metaheuristic algorithms has gained attention for applications where the resistivity distribution exhibits sharp boundaries or extreme contrasts that challenge conventional smoothness-constrained inversions. PSO-based 2D ERT inversion has been demonstrated for hydrogeological applications, with results showing improved delineation of resistivity boundaries compared to Occam-style inversions [54]. GA-based approaches have been applied to characterize waste disposal sites and assess groundwater contamination plumes.

GPR: FWI of GPR data using metaheuristic algorithms has been explored for near-surface characterization. DE-based GPR inversion has shown promise for estimating permittivity and conductivity profiles in the shallow subsurface, with applications in archaeology, engineering, and environmental site investigation [55]. PSO has been applied to optimize GPR survey parameters and interpret radargram attributes for infrastructure assessment.

CSEM: marine CSEM inversion for hydrocarbon exploration has been addressed using SA and GA approaches. The high computational cost of 3D CSEM forward modeling (typically 30–120 minutes per model evaluation) severely limits the number of forward simulations feasible within a metaheuristic framework, necessitating either dimensionality reduction or surrogate model approximations. Recent hybrid approaches combining PSO with local Gauss-Newton optimization have shown promising results for 1D and 2D CSEM inversion [56].

Seismological applications of metaheuristic algorithms extend beyond exploration seismology to encompass earthquake seismology, where the problems of hypocenter location, focal mechanism determination, seismic hazard analysis, and earth structure estimation present unique optimization challenges. The review identified 22 publications (11.8%) in this category.

Earthquake hypocenter location: the determination of earthquake hypocenter coordinates (latitude, longitude, depth) and origin time from arrival-time data is a fundamental seismological inverse problem. While linearized methods (e.g., Geiger's method) work well for events recorded by dense networks with good azimuthal coverage, they can fail for poorly constrained events, particularly regarding depth estimation. SA was applied to earthquake location by Billings et al. [57], demonstrating robust performance for events with limited station coverage. PSO-based earthquake location algorithms have shown improvements over grid-search methods, particularly for events in complex velocity structures where the travel-time residual surface contains multiple minima. The Oct-Tree method of Lomax et al. [58] and the NA approach of Sambridge and Kennett [59] represent related stochastic approaches that have been widely adopted in seismological practice.

Focal Mechanism Determination: the determination of earthquake focal mechanisms (strike, dip, and rake angles) from seismic waveform data or first-motion polarities has been addressed using GA and SA approaches. The nonlinear relationship between source parameters and observed waveforms, combined with the presence of multiple acceptable solutions due to data limitations, makes this problem naturally suited to global optimization. Snoke [60] applied GA to focal mechanism determination, while Sokos and Zahradník [61] demonstrated SA-based moment tensor inversion. Recent applications have employed DE for regional moment tensor inversion with improved uncertainty quantification through ensemble analysis of near-optimal solutions.

Seismic Hazard Analysis: metaheuristic algorithms have been applied to optimize ground motion prediction equations (GMPEs), calibrate seismicity rate models, and optimize seismic network configurations. Kayhan et al. [62] applied Harmony Search to the optimization of structural design under seismic loading, while GA-based approaches have been used to optimize attenuation relationship coefficients. Sun and Qi [63] recently proposed a Generalized Chaotic PSO framework for Ground Motion Inversion (GCPSO-GIT), achieving a site effect coefficient of variation of only 12% compared to 35% for standard PSO-GIT, with a median source stress drop of approximately 42 bar and peak ground acceleration predictions within 3.6% of observed values.

The application of metaheuristic algorithms to well log analysis and petrophysical parameter estimation represents a growing area with significant practical implications for reservoir characterization. The review identified 18 publications (9.7%) in this domain, addressing problems including mineral volume estimation, porosity-permeability prediction, fluid saturation estimation, rock physics model calibration, and electrofacies classification.

Mineral volume estimation: the inversion of well log data (gamma ray, density, neutron porosity, photoelectric factor, sonic) for mineral composition and porosity is a constrained nonlinear optimization problem. The constraint that mineral fractions must sum to unity and remain non-negative naturally lends itself to bound-constrained metaheuristic optimization. PSO has been applied to estimate mineral volumes from multi-log data, outperforming conventional linear programming approaches when nonlinear mixing rules are employed [64]. DE-based mineral inversion has demonstrated improved robustness to log measurement errors and the ability to incorporate complex mineral-fluid interaction models.

Rock physics inversion: the estimation of rock physics properties (e.g., Hertz-Mindlin parameters, cementation exponents, critical porosity) from seismic and well log data involves calibrating theoretical rock physics models to observed data. GA and SA have been applied to optimize the parameters of Gassmann's fluid substitution equations and granular media theories, enabling prediction of elastic properties under varying saturation and pressure conditions [65]. PSO-based rock physics inversion has shown promise for predicting reservoir quality from seismic attributes, with applications in carbonate reservoir characterization. Facies Classification: metaheuristic algorithms have been employed to optimize the parameters of clustering and classification algorithms for electrofacies determination from well logs. GA-optimized neural network architectures and PSO-tuned support vector machine hyperparameters have shown improvements in facies prediction accuracy compared to default parameter settings [66].

The optimization of geophysical survey design including sensor placement, source-receiver geometry, acquisition parameters, and survey cost minimization represents a combinatorial optimization problem naturally suited to metaheuristic approaches. Although fewer in number (9 publications, 4.8%), these applications address a practically important challenge: maximizing the information content of geophysical measurements while minimizing acquisition costs and logistical complexity.

Optimal sensor placement: the placement of seismometers, gravity stations, magnetometers, or electrodes to maximize the spatial resolution of the resulting geophysical images is a combinatorial optimization problem. ACO has been applied to optimize seismic receiver placement for microseismic monitoring of hydraulic fracturing, maximizing the expected number of detectable events while minimizing network cost [67]. GA-based survey design has been demonstrated for gravity monitoring networks, where the objective is to maximize the sensitivity to anticipated subsurface mass changes (e.g., CO₂ plume migration or groundwater level fluctuations).

Acquisition geometry optimization: the design of seismic acquisition geometries including source and receiver line spacing, offset distribution, and azimuthal coverage has been addressed using PSO and GA. These approaches optimize proxy measures of acquisition quality, such as fold distribution uniformity, offset-azimuth coverage, and spatial sampling completeness, subject to operational constraints including access restrictions, cost limits, and equipment availability. DE-based optimization of marine seismic streamer configurations has shown improvements in subsurface illumination for complex geological targets.

Survey cost minimization: multi-objective metaheuristic approaches, particularly NSGA-II (non-dominated sorting GA II), have been applied to optimize the trade-off between survey information content and acquisition cost. These approaches generate Pareto-optimal survey designs that allow geophysicists to select the most cost-effective configuration for their specific objectives and budget constraints [68].

A central objective of this review is to provide a quantitative cross-domain comparison of metaheuristic algorithm performance for geophysical inverse problems. Drawing from the 186 reviewed publications and normalizing performance metrics to enable comparison across different geophysical modalities, Table 8 summarizes the best-performing algorithms for each major inverse problem type.

The analysis reveals several consistent patterns. First, PSO emerges as the most versatile algorithm, achieving top-two rankings in five of the eight problem categories examined. This is attributed to its simple implementation, rapid convergence for moderate-dimensional problems (M<50), and effective balance between exploration and exploitation. Second, DE demonstrates consistently strong performance, particularly for problems with continuous parameter spaces and moderate-to-high dimensionality, outperforming PSO in gravity basement estimation and pre-stack seismic inversion. Third, SA (and its variants VFSA and ASA) maintains advantages for high-dimensional problems where population-based methods become computationally prohibitive, as the single-solution trajectory of SA requires only one forward model evaluation per iteration. Fourth, newer algorithms (GWO, WOA, HHO) show competitive performance in specific applications but lack the extensive validation and parameter sensitivity analysis available for established methods.

Note: NS=Not statistically significant. Statistical significance determined by Wilcoxon rank-sum test across multiple independent runs. Misfit improvement calculated relative to the third-best performing algorithm in each category.

One of the most significant trends identified in this review is the increasing adoption of hybrid metaheuristic approaches that combine the global exploration capability of stochastic algorithms with the rapid local convergence of gradient-based methods. These hybrid strategies address the fundamental limitation of standalone metaheuristics: while they excel at locating the basin of attraction of the global optimum, they are often inefficient at fine-tuning the solution within that basin, requiring many additional forward model evaluations that would be unnecessary with local optimization. The most common hybridization strategy involves running a metaheuristic algorithm for a fixed number of iterations (the global phase) to identify promising regions of the model space, followed by local gradient-based refinement (the local phase) from the best solution found. The LM algorithm is the most frequently used local optimizer in these hybrids due to its robustness and adaptive damping capabilities. Alternative local methods include the Gauss-Newton algorithm, conjugate gradient methods, and the Occam inversion framework. A second hybridization strategy involves embedding local search within the metaheuristic iteration loop, where selected individuals in the population undergo local optimization at regular intervals (memetic algorithms or Lamarckian evolution).

Monte Carlo methods, particularly MCMC and the NA, occupy a distinct niche in geophysical inversion as they provide rigorous posterior probability distributions rather than point estimates. The comparison between metaheuristic optimization and Monte Carlo sampling methods is nuanced, as they serve complementary purposes: metaheuristics seek the optimal model (maximum likelihood or maximum a posteriori estimate), while MCMC methods characterize the full posterior distribution and provide formal uncertainty quantification.

Sambridge [33], [45] introduced the NA as a hybrid approach that bridges optimization and Bayesian inference. The NA first performs a stochastic search (NA-search) to identify high-likelihood regions of the model space, then resamples the ensemble using importance weights (NA-resampling) to approximate the posterior distribution. Comparisons in the reviewed literature suggest that for point estimation (finding the best-fit model), PSO and DE typically converge 2–4 times faster than MCMC while achieving comparable or superior misfit values. However, MCMC provides statistically rigorous uncertainty estimates that are lacking in standard metaheuristic implementations. Recent work on PSO-based uncertainty quantification, where the ensemble of personal best solutions is used to approximate the posterior, shows promise but lacks the formal convergence guarantees of MCMC [71].

To provide a rigorous statistical comparison across the reviewed studies, a Friedman ranking test was performed on algorithm performance across the six geophysical application domains. Each algorithm was ranked within each domain based on a composite performance metric combining normalized misfit, convergence speed, and consistency (coefficient of variation across multiple runs)O F = 42.8 (p < 0.001) confirmed that significant differences exist among algorithms across geophysical domains. Nemenyi post- O R formed a statistically indistinguishable top tier, significantly outperforming all other algorithms, 2) VFSA and GA formed a second tier, not significantly different from each other but significantly worse than PSO and DE, 3) GWO and ABC formed a third tier with competitive but variable performance and 4) WOA, HHO, MFO, BA, SCA, and ACO showed the most variable performance, with effectiveness strongly dependent on the specific geophysical problem structure. These rankings should be interpreted cautiously, as they aggregate across diverse problem types, dimensionalities, and study quality levels. Fig. 2 presents the Critical Difference (CD) diagram derived from the Friedman ranking test followed by the Nemenyi post-hoc statistical analysis for twelve metaheuristic algorithms evaluated across six geophysical application domains. The horizontal axis represents the average Friedman rank, where lower values indicate superior overall performance. PSO achieved the best average rank (1.83), closely followed by DE (2.17), and the absence of statistically significant differences between these two algorithms is indicated by the connecting horizontal bar. VFSA (3.50) and GA (3.83) form a second statistically equivalent performance group, demonstrating competitive but moderately lower effectiveness compared to PSO and DE. GWO (5.17) and ABC (5.50) constitute a third performance cluster with acceptable but more variable behavior across problem domains. The remaining algorithms, including WOA, HHO, MFO, BA, SCA, and ACO, occupy lower-ranking positions between 7.0 and 10.5 and exhibit overlapping confidence intervals, indicating inconsistent performance depending on problem characteristics and dimensionality. The CD O ORO statistical analysis confirms that PSO and DE represent the most consistently reliable general-purpose optimization algorithms for geophysical inversion problems among the methods reviewed in this study.