Introduction
Induced polarization tomography (IPT) is a non-invasive geophysical method for subsurface imaging based on changes in the polarization parameter, which is used in various fields such as subsurface resource exploration and environmental studies. The main goal of induced polarization tomography data inversion is to infer subsurface physical properties from noisy data, which is often accompanied by various uncertainties. These uncertainties are divided into two main categories, including inherent data noise and incomplete knowledge of the physics of the problem. In this study, the uncertainty of polarization models resulting from induced polarization tomography data inversion using moving block bootstrap (MBB) resampling methods is considered. This method helps to analyze uncertainty in models by generating a set of models with a defined level of fit to the observed data. The moving block bootstrap is specifically designed to improve sampling accuracy in data where there is a possibility of correlation and dependence between data samples. This approach is used to assess the uncertainty of model parameters under different conditions, both in terms of synthetic data and field data. The results of this research can help improve geophysical analyses and decision-making related to predicting subsurface conditions. Numerical modeling on both synthetic and field data shows that increasing the complexity of the subsurface model and decreasing the model’s sensitivity to the data leads to greater uncertainty in the inverted model. Additionally, the type of regularization term and the choice of regularization parameters were found to influence uncertainty quantification. In this study, we developed a software package in the MATLAB programming environment that includes comprehensive modules for forward modeling, effective resampling using the aforementioned approach, and inverse modeling. The software has been tested on various synthetic examples and field data. Our numerical tests involve two synthetic models with specified chargeability and resistivity values, designed for reconstructing areas prone to mineral deposits and valuable metal and sulfide minerals covered by alluvial overburden. For real-world examples, two profiles from two different regions in northwestern and southeastern Iran were selected.The first study area includes a 2000 m geoelectrical profile, and the second area, with an 800 m geoelectrical profile and the presence of exploratory boreholes, provides robust evidence for evaluating and testing the proposed approach.
Materials and Methods
This research follows a step-by-step analytical workflow including: forward modeling of induced polarization response, generation of synthetic data with controlled noise, inversion of resistivity and induced polarization data, implementation of the moving block bootstrap resampling algorithm, and finally statistical analysis of results.Forward Modeling: To simulate the earth's response to time-domain induced polarization, the chargeability perturbation model proposed by Seigel (1959) was used. Apparent chargeability (η_a) is calculated from:
Eq. 1)


where σ is conductivity, η is chargeability, V_p is primary voltage, V_s is secondary voltage, and F_dc is the DC resistivity forward modeling operator. This equation reveals that the apparent chargeability is computed based on two consecutive applications of the resistivity forward operator for two conductivity models, that is σ and σ(1-η).
Inversion: For data inversion, the non-linear approach of Oldenburg and Li (1994) was used. The objective function, consisting of the weighted sum of data misfit (Φ_d) and model stability (Φ_η), is defined as:
Eq. 2)
.

Here, η is the chargeability model, α is the damping factor, W_d is the data weighting matrix (inverse of data error), W_m is the model smoothing matrix (combination of horizontal and vertical derivatives), and η_apr is the prior chargeability model. The inversion starts with an initial model η^0 obtained from the geometric mean of apparent chargeability data. The IP sensitivity matrix is updated at each iteration, and the iterative numerical method continues until convergence. Moving Block Bootstrap (MBB) Resampling: To quantify uncertainty, the MBB resampling algorithm as a non-parametric method is used. Unlike classical bootstrap that assumes data independence, this method preserves the correlation and dependence structure of the data points. The bootstrap resampling algorithm has found numerous practical applications in the mathematical literature for statistical inference, such as estimating distribution functions, conducting hypothesis tests, and identifying models. Despite its proven effectiveness, it has been seldom applied for quantitatively assessing uncertainty intervals in geophysical models. Following data resampling, the inversion process is performed independently for each data realization (bootstrapped data). Finally, among the inverted models, those with an appropriate fit level (i.e., χ^2≤m+√2m, where m is the number of data) are selected as equivalent models, and their mean and standard deviation are calculated. The mean represents the stable model structure, while the standard deviation represents the uncertainty map of the inverted chrgeability models.
Results and Discussion
Results from synthetic modeling and field data sets indicate that the Moving Block Bootstrap (MBB) resampling approach is an effective tool for quantifying uncertainty in the inversion of induced polarization tomography data. In both designed synthetic models, it was observed that the mean model derived from 1000 resampling runs reconstructed the main subsurface structures (including the location of anomalies, layer boundaries, and alteration zones) with reasonable accuracy compared to a single inversion model. The standard deviation maps from the uncertainty analysis systematically identified the areas with the highest model variability. In the first synthetic model, the highest uncertainty was observed at the boundary between the surface layer and the host medium, as well as at the edges of the anomaly. In the second, more complex synthetic model (which includes a fault structure and three zones with different properties), uncertainty increased not only at the overburden-fault boundary but also within the environment exhibiting higher chargeability. These patterns indicate that uncertainty is mainly concentrated in areas with sharp gradients in physical parameters, as well as in parts of the model where data sensitivity to the parameters decreases. This finding is consistent with the theoretical foundations of the ill-posed inverse problem and demonstrates that the MBB method can naturally highlight areas requiring supplementary surveys or more detailed analyses. In the field data, a significant agreement was observed between high-chargeability zones in the inverted models and the mineralization intervals recorded in exploratory boreholes (N23 and N12). This depth and lateral consistency confirms the validity of the proposed method under real conditions. Particularly in Profile 2 from the Khusf area, the coincidence of zones with chargeability exceeding 14-15 mV/V and resistivity below 20 Ωm with the sulfide mineralization logs demonstrates the high efficiency of the
method in distinguishing alteration and mineralization zones from barren rocks. It was also observed that the highest uncertainty lies in the transition zones between resistive rocks and altered zones, which matches the lithological changes reported in the boreholes. Importantly, the MBB method not only helps identify mineralization zones but also reveals the boundaries of high uncertainty, which represent geologically high-risk areas. In the first field dataset (Meshgin-Shahr area), which lacked exploratory boreholes, the uncertainty analysis showed that the western part of the survey profile had the highest standard deviation. This section had previously shown different behavior in pseudo-sections, characterized by low resistivity and high chargeability. The alignment of this high uncertainty with a possible alteration zone (argillic) suggests that even in the absence of direct control data, the standard deviation map can serve as an indicator for identifying complex areas requiring supplementary investigation. On the other hand, the results showed that the type of regularization term and its coefficients have a direct impact on the estimated uncertainty. The use of the L_2 norm together with horizontal and vertical smoothing matrices led to more continuous models and a relative reduction of uncertainty in homogeneous areas, but at the same time somewhat smoothed the sharp boundaries of geological structures. These observations underscore the necessity of optimizing regularization parameters according to the study objective (distinguishing continuous versus discrete structures). In summary, the results of this discussion confirm that the MBB resampling approach, beyond merely quantifying uncertainty, provides the following capabilities: identifying high-sensitivity areas requiring further investigation, indirect validation of inversion results in the absence of control data, and providing a basis for decision-making in the design of supplementary surveys. These features are of high practical value, especially in mineral exploration and environmental studies where the cost of drilling and additional surveys is very high.
Conclusion
In this study, a comprehensive framework for quantifying uncertainty in induced polarization tomography data inversion using the non-parametric Moving Block Bootstrap (MBB) resampling method was presented and implemented. The developed algorithm in MATLAB is capable of forward modeling, controlled noise application, 2D inversion of resistivity and chargeability, and generation of multiple resampled datasets with adjustable block length, removal percentage, and number of iterations (1000 times). Results from synthetic models showed that standard deviation maps clearly identify high-uncertainty areas, which mainly correspond to geological boundaries, anomaly edges, and zones with sharp gradients in physical parameters. For field data, significant agreement between high chargeability zones in inverted models and mineralization intervals recorded in exploratory boreholes (Khusf area) was observed, confirming the validity of the proposed method. Furthermore, in the Meshgin-Shahr area (lacking boreholes), the uncertainty map successfully identified the western part of the profile as a complex area requiring supplementary investigation.