Iterative geostatistical electrical resistivity tomography inversion

Electrical resistivity tomography (ERT) is a geophysical method used to create an image of the subsurface due to its sensitivity to porosity, water saturation, and pore fluid salinity. This geophysical method has been widely applied in the investigation of mineral and groundwater resources, as well as in archaeological, environmental, and engineering studies. The prediction of subsurface properties, such as electrical conductivity, from measured ERT data requires solving a challenging geophysical inversion problem. This work proposes an iterative geostatistical resistivity inversion method using stochastic sequential simulation and co-simulation as stochastic model perturbation and update techniques. Electrical resistivity models are generated conditioned to a target histogram, often retrieved from available resistivity borehole data, and assuming a spatial continuity pattern described by a variogram model. From the electrical resistivity models, a finite-volume approximation of Poisson’s equation is used to compute synthetic ERT data. The misfit between predicted and observed data drives the convergence of an iterative procedure and conditions the co-simulation of new models in the subsequent iterations. This methodology is applied to a two-dimensional synthetic case, and a set of two-dimensional profiles obtained from an ERT survey carried out in southern Portugal. In both application examples, the final models predict ERT data that match the observed ones while reproducing borehole data and imposed variogram models. The results obtained in both data sets are compared against a commercial deterministic ERT inversion methodology, showing the ability of the proposed method to model small-scale variability and assess spatial uncertainty.


Introduction
Electrical resistivity tomography (ERT) is a geophysical method used to predict the spatial distribution of subsurface electrical resistivity (e.g., Parasnis 1986;Telford et al. 1990;Reynolds 2011).Since electrical resistivity is directly related to rock type, porosity, ionic strength of the pore fluids, and surface conductivity of geologic materials (Sumner 1976;Sharma 1997), ERT is widely used in hydrological studies (Page 1968;Wilson et al. 2006), mineral exploration (Bauman 2005;White et al. 2001), archaeological prospection (Griffiths and Barker 1994;Tsokas et al. 2008), and environmental and engineering studies (Chambers et al. 2006;Rucker et al. 2010).ERT data are collected by establishing an electrical potential difference between two current electrodes.An electrical current is injected into the ground, and the resulting potential distribution is measured at many pairs of potential electrodes (Griffiths and Barker 1993).The observed measurements are then converted into apparent resistivity, which represents a weighted average of the resistance of earth materials to current flow, providing a smooth representation of the true subsurface spatial distribution of electrical resistivity (Loke et al. 2013).Apparent resistivity enables a qualitative prediction of the electrical parameters Published in the special issue "Geostatistics and hydrogeology" of the subsurface, but it is not sufficient to predict and capture the true spatial distribution and variability of the subsurface electrical resistivity.Apparent resistivity models distort the real subsurface characteristics as these data correspond to volumetric measurements highly dependent on the type and configuration of the acquisition (Dahlin and Zhou 2004;Saydam and Duckworth 1978).
To predict the true subsurface electrical resistivity spatial distribution, observed apparent resistivity needs to be inverted (Loke 2002).Due to measurement errors in the acquisition, noise contamination, and incomplete data coverage, ERT inversion is an ill-posed, nonlinear problem with a nonunique solution (e.g., Tarantola 2005).Multiple solutions imply uncertainty about the prediction obtained.Hence, an accurate assessment of model uncertainty is fundamental to properly interpreting the predictions and for well-informed decision-making.
The relationship between observed geophysical data and model parameters can be mathematically described as where m represents the vector model parameters to be predicted (i.e., electrical resistivity), d obs corresponds to the vector of observed data (i.e., apparent resistivity), F is the forward operator that maps the model into the data domain, and e is a vector that represents the discrepancies arising from measurement errors and the assumptions made during the data processing.
Classical ERT inversion methods are deterministic.A deterministic inversion searches for a single earth model (i.e., the expected model) able to predict ERT data with an acceptable fit to the observed data, satisfying any other imposed constraints such as being consistent with an initial model built from the a priori knowledge about the subsurface geology.In deterministic inversion methods, the solution minimizes an objective function consisting of a regularized weighted least squares formulation, in which the search is usually conducted using iterative gradient-based methods (e.g., Ellis and Oldenberg 1994;LaBrecque et al. 1996;Loke and Barker 1996;Pidlisecky et al. 2007).As deterministic inversion predicts a single solution, it does not provide insights into the degree of uncertainty associated with the inversion results.Several works have investigated the use of geostatistical priors to regularize and impose a given spatial continuity pattern to the predicted models (e.g., Hermans et al. 2012;Jordi et al. 2018;Bouchedda et al. 2017).Hermans et al. (2016) propose the prediction-focused approach (PFA) forecast the spatiotemporal change hydrogeological properties using electrical resistivity tomography.Linde et al. (2015) review the most common methods to include geological realism in hydrogeological and geophysical inverse modelling.
(1) d obs = F(m) + e Alternatively, stochastic geophysical inversion methods search for multiple subsurface models (e.g., of electrical conductivity) that predict ERT data that fit equally well with the observed ERT data.A variety of stochastic methods are described in the literature, but in general, they are divided into two main groups of techniques: Bayesian inversion and stochastic optimization algorithms (e.g., Tarantola 2005; Gloaguen et al. 2005;Giroux et al. 2007;de Pasquale et al. 2017de Pasquale et al. , 2019;;de Pasquale and Linde 2017).
In Bayesian inversion methods, a joint posterior probability distribution for all model parameters is used to describe the solution to the inverse problem.The posterior distribution is obtained using a likelihood function built on the available data sources, which updates a prior distribution for the model parameters.Zhang et al. (1995) proposed an inversion method to maximize model parameters' posterior probability density function.This method's implementation for ERT relies on assumptions about the spatial covariance of the resistivity parameters and Gaussian distributions for data errors and model parameters.Mosegaard and Tarantola (1995) described a statistical approach reformulated as a Bayesian inference problem using Markov chain Monte Carlo (McMC) and the Metropolis algorithm sampling method.In their work, the posterior distribution combines physical models and available prior information with new information obtained through direct measurement of the subsurface.
Alternatively, stochastic optimization algorithms methods to predict hydrogeological properties approximate the posterior distribution.A comparison between three global optimization methods is provided in Barboza et al. (2018).The most cited works include ERT inversion based on McMC to assess the posterior distribution of the model parameters shown in de Pasquale (2017), de Pasquale and Linde (2017), de Pasquale et al. (2019) and Aleardi et al. (2020).Chen and Zhang (2006) propose an ensemble Kalman filter for providing updated estimates of model parameters and model state variables, such as hydraulic conductivity and pressure head and their uncertainty.Arboleda-Zapata et al. (2022) propose a workflow to analyze ensembles of predicted electrical resistivity models.
In iterative geostatistical geophysical inversion methods (e.g., Azevedo and Soares 2017;Grana et al. 2021), the model parameters are considered as a realization of a random function.In this context, the model parameter space is perturbed and updated using stochastic sequential simulation and co-simulation coupled with a global optimizer.The optimization is driven by the mismatch between observed and predicted synthetic data.At the end of the iterative procedure, a set of subsurface models representing the posterior probability distribution are obtained.The uncertainty of the predicted models can be assessed, for example, by computing the pointwise inter-quantile range of the set of inverted models.The application of these methods to invert ERT data is still limited.
Dealing with ERT data, Yeh et al. (2002) describe a sequential geostatistical ERT inversion method that uses spatial covariance matrices to include prior knowledge about general geological structures.The method uses well-log data to constrain the solution and a successive linear estimator to find an optimal model.Feyen and Caers (2006) employed multiple-point geostatistics to characterize the hydrofacies architecture of complex geological settings, using a training image designed to reflect the prior geological knowledge.They also used a spatial covariance and a multi-Gaussian random function to model the intra-facies variability of the hydraulic properties.Mariethoz et al. (2009) used truncated pluri-Gaussian simulation to assess contaminant migration in highly heterogeneous porous media.Truncated pluri-Gaussian simulation attempts to create maps of categorical values by truncating at least two underlying multi-Gaussian simulations.Hörning et al. (2020) presented a geostatistical approach for the inversion of ERT data based on "random mixing"; in this technique, realizations of conductivity fields are constructed by combining random fields that have the spatial correlation of conductivity.Lochbühler et al. (2014) propose a method to condition the generation of subsurface models with multiple-point statistics with tomographic images.
This paper proposes an alternative iterative geostatistical resistivity inversion method based on direct sequential simulation and co-simulation (Deutsch and Journel 1998).The available resistivity borehole data are used to model the spatial continuity patterns of subsurface geology as given by a variogram model and to condition the generation of electrical resistivity models locally.When borehole data are not available, variogram models and the target distribution of electrical resistivity retrieved from analogues can be imposed.The similarity coefficient between observed and predicted data at a given iteration drives the convergence of the iterative procedure and the assimilation level of the observed ERT data from iteration to iteration.
The proposed geostatistical ERT inversion methodology is applied to two-dimensional (2D) synthetic and real data sets.The real application example consists of a set of 2D profiles obtained from an ERT survey carried out in the Alentejo region in Portugal, which was designed to model the groundwater system of the region.All predicted models generate synthetic ERT data similar to the observed ones, reproducing both the borehole data and the imposed variogram models.In addition, the predicted models show the ability of the proposed inversion method to characterize the spatial uncertainty of the model parameters.The results obtained in both application examples are compared against a conventional deterministic inversion methodology available in commercial software (RES2DINV) (Loke 2010).
The next section details the proposed methodology.This is followed by synthetic and real case applications, including a detailed description of the data sets.The results are discussed in the subsequent section before the main conclusions.

Methodology
The proposed iterative geostatistical ERT inversion method aims to predict the spatial distribution of subsurface electrical resistivity from recorded ERT data (i.e., apparent resistivity pseudo-sections).This method encompasses three main steps-model generation, generation of synthetic ERT data, and a stochastic update (Fig. 1).Each step is described in detail in the following.

Model generation
Electrical resistivity model parameters are generated using direct sequential simulation (DSS) during the first iteration and direct sequential co-simulation (co-DSS) in the following iterations (Soares 2001; Azevedo and Soares 2017).At each iteration, a set of Ns electrical resistivity models are generated, accounting for the spatial continuity given by a variogram model (i.e., a spatial covariance matrix), considering the local distribution that the attribute should have at each location and conditioned to the available resistivity borehole data.The variogram model and the local probability distribution can be estimated from the existing borehole data and/or inferred from expert knowledge.
Briefly, in direct sequential simulation, each location of the simulation grid is visited sequentially following a random path.At each visited location, a value of the original variable is drawn from a probability distribution function based on a simple kriging estimate using observed data (i.e., direct observations) and previously simulated values within a given neighborhood (Deutsch and Journel 1998;Soares 2001).The simple kriging estimate and variance are used to build an auxiliary probability distribution function from the distribution of the experimental data set.The stochastic sequential simulation finishes after all the locations of the random path are visited.Each time the simulation runs, an alternative model is generated (i.e., a geostatistical realization) as the random path changes, and so they change the previously simulated data at each location.

Generating synthetic ERT data
The forward model implemented in the iterative geostatistical ERT inversion method described herein was developed by Pidlisecky and Knight (2008).During the iterative procedure, this 2D forward model is used to compute Ns synthetic apparent resistivity models from the previously generated Ns electrical resistivity geostatistical realizations.Using a 2D forward model might represent a limitation when computing the ERT response from highly complex geological settings, as the injected electrical current into the subsurface flows three-dimensionally through preferential paths that could circumvent resistive structures present in a 2D representation.In these cases, alternative 3D forward models could be used, but the computational costs of the proposed methodology would increase.A summary of the main principles followed by Pidlisecky and Knight (2008) is provided in the following.
In ERT surveys, a series of voltage measurements are obtained in response to a series of known input currents.Poisson's equation can be used to describe the electric potential field generated when a current passes across an electrode dipole where σ is the electrical conductivity , δ is the Dirac delta function, and r + and r -are the locations of the positive and negative current electrodes, respectively.To numerically solve Eq. ( 2) for the electric potential, ϕ, numerical gradient, and divergence approximations are required.Following Pidlisecky and Knight (2008), once numerical finite difference operators have been derived for gradient and divergence, Eq. ( 2) can be written in matrix notation as where D is the divergence matrix, S(σ) is a diagonal matrix containing the electrical conductivity values, G is the gradient matrix, φ is a vector of electric potentials, A(σ) is the combined forward operator, and q is a vector containing the current electrode pairs.Equation ( 3) is solved to yield the potential field Equation (4) results in a vector of electric potential values for the cells in the model.Knowing the survey potential electrode locations, potential differences can be calculated across each measurement pair.These measurements are then multiplied by the geometric factor K to provide apparent resistivities The geometric factor (K) depends on the arrangement of the four electrodes (i.e., depends on the distance between each electrode and the measurement).K is given by where r C1-P1 is the distance between current electrode C1 and potential electrode P1, r C1-P2 is the distance between current electrode C1 and potential electrode P2, r C2-P1 is the distance between current electrode C2 and potential

Stochastic update
At the end of each iteration, the stochastic update of the electrical resistivity models is performed using a data selection procedure that controls the assimilation degree of the observed ERT data.For a given iteration, and for the set of Ns simulated electrical resistivity models, Ns synthetic apparent resistivity models are computed using the forward model described previously.The computed apparent resistivities are locally compared against the observed one in terms of a similarity coefficient (S) using a nonoverlapping moving window that visits all the inversion grid locations where x and y are the observed and synthetic apparent resistivity, respectively.N is the number of observations used in the calculations.The moving window does not need to be square, and its width and height are randomly generated at the beginning of each iteration to avoid biasing the results from iteration to iteration.Alternative similarity coefficients could be used as long as they are bounded between -1 and 1 with a similar meaning to Pearson's correlation coefficient.This assumption is required due to the use of S in the stochastic sequential co-simulation of new models in the subsequent iteration.
At each moving window location, the samples of electrical resistivity corresponding to a given geostatistical realization from which the largest similarity coefficient originated are stored together with the similarity coefficient in two auxiliary arrays, which are used as conditioning information in the subsequent iteration.
In the new iteration, a new set of Ns models is co-simulated using both auxiliary arrays as secondary variables.The magnitude of S determines the variability of the new ensemble of electrical resistivity co-simulated models.The higher the similarity coefficient, the less variable the ensemble will be (i.e., the higher the assimilation of the observed apparent resistivity data).S is similar to Pearson's correlation coefficient but is sensitive to the amplitude mismatch between signals.The iterative process finishes when the similarity coefficient, computed over the entire domain, is above a given threshold, or if a given number of predefined iterations is reached.
During the entire iterative procedure, each electrical resistivity model generated with DSS and co-DSS reproduces the observed data at their locations, the probability distribution function of electrical resistivity, and the variogram model imposed during the stochastic sequential (co-)simulation.
The variogram model adopted for the inversion depends on the data availability and will condition the geological plausibility of the predicted subsurface models.The proposed iterative geostatistical ERT inversion method can be summarized by the following sequence of steps (illustrated in Fig. 1

Application examples
The proposed iterative geostatistical ERT inversion methodology was applied to 2D synthetic and real data sets.The synthetic application example acts as proof of concept of the proposed methodology and is compared against a commercial deterministic inversion solution.The results obtained with the real application example that considers realistic noise levels are compared against a conventional deterministic inversion methodology to analyze its advantages and disadvantages.

Synthetic application example
The synthetic application example shown herein is inspired by a laboratory experiment conducted in a sandbox by Citarella et al. (2015) and Chen et al. (2018).
In this experiment, a laboratory sandbox is filled with homogeneous porous material (i.e., glass beads) and has an impermeable barrier positioned in the middle top of the sandbox.Within the sandbox, pollutant dispersion in a groundwater system is simulated by injecting a tracer solution into the porous medium and controlling the head level.A photometric method is used to monitor the plume evolution in time.This experiment mimics a typical groundwater system recharged by natural rainfall entering the soil profile and leaching into deeper soil layers.Due to intensive agricultural or industrial activities, the leachate leaving the soil profile and entering the aquifer may contain concentrations of toxic substances.Once these substances have entered the aquifer, they can be transported over large horizontal distances, thus contaminating large parts of the aquifer.In the case of groundwater contamination, it is important to understand how the toxic substances are dispersing so that proper mitigation actions can be taken.The inversion results illustrated herein aim at assessing the potential of the proposed inversion method to detect contamination plumes.
The electrical resistivity reference data set used in this work represents a snapshot of the system described previously with the plume already dispersed under the impermeable barrier (Fig. 2a).Plume spread is apparent by the V-shaped low resistivity feature in a high resistive background (i.e., the glass beads filled with freshwater).The vertical impermeable barrier induces the V-shape.The sandbox is 90 cm long by 18.2 cm high.Tracer dispersion occurs from right to left.The impermeable barrier is observed as a vertical low resistivity feature starting from the top of the model until a depth of about 5.6 cm and positioned at a horizontal distance of 43.5 cm from the left border.The 2D inversion grid consists of 60 × 1 × 13 cells for the i-, j-, and k-directions, respectively.
The reference apparent resistivity (Fig. 2b) was numerically computed considering a Wenner-Schlumberger acquisition array (e.g., Loke 2002) composed of a total of 31 electrodes spaced every 3 cm and solving the forward model to yield the potential field following Pidlisecky and Knight (2008).This apparent resistivity was used as true geophysical data during the application of the proposed methodology.The same forward model used to calculate the true apparent resistivity field was used as part of the inversion.This approach assumes that no uncertainty is considered in the forward model, which might be a strong assumption in real case applications with complex geology settings.
To apply the proposed geostatistical resistivity inversion to the reference data set, the true electrical resistivity field was sampled at two boreholes on both sides of the impermeable barrier.The position of the boreholes can be seen in Fig. 2. The borehole data were used as experimental data to condition the generation of models during the iterative geostatistical inversion.As two boreholes are being considered as conditioning data, the spatial continuity pattern of both horizontal and vertical directions was estimated directly from the true electrical resistivity, as represented by a 2D global variogram model (Fig. 3).In this synthetic application example, it is assumed that there is no uncertainty in the spatial continuity pattern imposed during the iterative procedure (i.e., no uncertainty on the variogram model).Also, to reduce the complexity of the synthetic data set, the simulation and inversion area is limited to the region where the apparent resistivity exists (Fig. 2b).
The experimental variograms in the horizontal and vertical directions were fitted with a spherical variogram model.The ranges used were 30 cm for the horizontal direction and 5 cm for the vertical one.
The iterative geostatistical inversion ran with 32 models of resistivity and for 6 iterations.These values were set after trial-and-error to make sure the iterative procedure converged.The evolution of the global S between reference and synthetic apparent resistivity is shown in Fig. 4. The models generated during the first iteration of the inversion procedure are characterized by a global S higher than 0.90.This means that the inversion problem is well characterized since a good convergence is reached at the early stages of the iterative procedure.This effect might be due to the relatively small size of the inversion grid versus the number of experimental data; also, the imposed variogram model is close to the true one.After the second iteration, the global S is higher than 0.95, reaching almost 1 at the end of the six iterations.As the stopping criterion for the iterative procedure, a fixed number of iterations was chosen, which was set after trial-and-error over a small portion of the area of interest.
Figure 5 shows the results of the inversion.In Fig. 5a, the reference apparent resistivity computed from the reference map in Fig. 2a is shown, whereas Fig. 5b shows the apparent resistivity from the realization that has the best fit (i.e., with the highest similarity coefficient), and Fig. 5c shows the apparent resistivity computed from the pointwise mean of all realizations generated during the last iteration of the iterative procedure.Next, Figs.5e,f show the similarity coefficients computed between the reference apparent electrical resistivity and that obtained from the best-fit electrical resistivity realization and the pointwise mean.There are small-scale differences around the areas where the tracer is being injected, the impermeable barrier is located, and at the plume front (black arrows in Fig. 5).These areas are characterized by high and abrupt resistivity contrasts, which have an impact on the quality of the inverted pseudo-sections.This effect is also noticeable in the local S values.
The best-fit electrical resistivity model (i.e., the one that produces apparent resistivity with the highest global S) and the pointwise mean of the electrical resistivity models predicted in the last iteration of the inversion process are shown respectively in parts b and c of Fig. 6.These models reproduce the overall spatial distribution of the pollutant dispersion observed in the reference model (Fig. 6a), but small differences are identified.Neither the plume V-shape nor the impermeable barrier are accurately reproduced, which is consistent with the small-scale differences identified between synthetic and reference apparent resistivity pseudo-sections in Fig. 5.It is a challenge for geostatistical simulation based on two-point statistics to reproduce small features such as the impermeable barrier or the curved shape seen in the plume spatial distribution.Alternative geostatistical methods such as multiple-point geostatistical simulation could perform better.Each electric resistivity model generated during the iterative procedure reproduces the histogram of the true model as retrieved from the borehole.This is an intrinsic property of direct simulation and co-simulation algorithms and is of great importance to ensure subsurface geological consistency.Figure 7 compares the reference histogram and the best-fit model histogram.
To assess the performance of the proposed geostatistical ERT inversion method, it is necessary to compare the results obtained with those from a deterministic inversion (Fig. 8).The deterministic inversion was obtained using RES2DINV (Loke 2010) with a default parameterization.The predicted pseudo-sections of apparent resistivity obtained at the last iteration of the deterministic inversion (Fig. 8a) can reproduce the main patterns of the true data.The local similarity coefficients between these data are high for the entire inversion grid (Fig. 8b).The predicted electrical resistivity (Fig. 8c) does reproduce the main V-shape of the electrical resistivity anomaly but is smooth and has a lower spatial resolution than the predicted model from the geostatistical inversion (Fig. 6).The small barrier at the top of the model is almost undetected by the predicted electrical resistivity model.Additionally, one of the advantages of using iterative geostatistical geophysical inversion methods is the ability to assess the spatial uncertainty associated with model predictions.Figure 8 shows the pointwise variance of electrical resistivity computed from the ensemble of models generated during the first and sixth iterations of the inversion procedure (Fig. 9a,b, respectively).It was assumed that the spatial distribution of electrical resistivity is only variable in the area with geophysical data (area inside the grey lines in Fig. 9).The remaining areas correspond to the constant high resistive background, so there is no variability.During the first iteration, the spatial uncertainty in the area of interest is only dependent on the location of the borehole data since no geophysical data have been assimilated yet.As expected, the variance increases with the distance from the experimental data.In the last iteration of the inversion process, the spatial uncertainty decreases drastically as the observed geophysical data are assimilated during the iterative procedure revealing areas where the match between observed and predicted data is less good (i.e., the predictions at these locations are more uncertain).

Real application example
The proposed iterative geostatistical ERT inversion methodology was applied to four 2D profiles obtained from an ERT survey carried out at the Neves-Corvo mining site (Alentejo region, Portugal).The survey aimed to characterize the spatial distribution of a groundwater system within mining premises.The full data set consists of a total of 22 apparent resistivity profiles.The application of the proposed geostatistical inversion is shown for four profiles that intersect each other, allowing for the assessment of the spatial coherency between predictions as each profile is inverted individually (Fig. 10).The predicted models with the proposed inversion methodology were compared against models inverted with a deterministic inversion method.
The ERT survey was performed with a Wenner-Schlumberger acquisition array configuration (e.g., Everett 2013).Table 1 summarizes the survey setup for the acquisition of each profile.
The histograms necessary for the stochastic sequential simulation and co-simulation were derived from previous ERT deterministic inversions due to the lack of wells drilled along the profile cross-sections.Therefore, the geostatistical simulation and co-simulation were not locally conditioned by any borehole information.Given the lack of direct observations and their spatial sparseness, the horizontal variogram models were retrieved from the sections of electrical resistivity obtained with a deterministic inversion approach provided by the data owner and adjusted for the expected geological knowledge of the area.This approach is similar to the workflow used in iterative geostatistical seismic inversion methodologies, where the horizontal variogram models are computed directly from the seismic data instead of the borehole data (Azevedo and Soares 2017) and were also proposed by Hermans et al. (2012).These approaches tend to overestimate the variogram ranges adding uncertainty to the imposed variogram model and it might result in unplausible predicted models.An alternative, commonly used in geostatistical seismic inversion, but not applied in this application example, is to optimize the variogram model by running small inversion pilot regions (i.e., mini-inversions).In these small inversion pilot regions, the inversion grid is divided into smaller regions, where multiple inversions run with different parameterizations.Then, the results are interpreted based on the geological knowledge of the study area and the parameters with the best results are used to invert the full inversion grid.The number of samples along the borehole path in the vertical direction is also limited and not able to capture the expected variability of the subsurface electrical conductivity.The limitations in estimating the spatial continuity pattern do impact the inverse solution.The resulting variogram models are shown in Table 2 and Fig. 11.The evolution of the global S between observed and synthetic apparent resistivity for the four inverted profiles is shown in Fig. 12.At the end of the inversion, the models generated for the different profiles reach a global S higher than 0.9.The models predicted during the first iteration produce synthetic ERT data with similarity coefficients between 0.6 and 0.85.The high convergence at an early   stage of the iterative procedure indicates that the electrical conductivity models generated in the first iteration, when there is no assimilation of the observed ERT data, might resemble the true subsurface geology.In cases where the variogram model is not geologically plausible, these global S values would be smaller.The predicted ERT data with the deterministic solution reached a S value above 0.95 for the four sections.The synthetic apparent resistivity data computed from the best-fit electrical resistivity model for profiles P15, P17, P19, and P20 are shown in Fig. 13.These apparent electrical resistivity pseudo-sections reproduce the main spatial patterns seen on the field apparent resistivity.The differences between synthetic and observed data are mainly identified at depth and in areas characterized by pronounced irregular shapes with abrupt apparent resistivity contrasts (black arrows in Fig. 13).The local S computed between observed and best-fit apparent resistivity for profiles P15, P17, P19 and P20 (shown in Fig. 14) confirms the lower quality of the inverted results in these areas.The predictions obtained for these locations are therefore uncertain as reflected by the pointwise variance computed from the ensemble of resistivity models computed during the last iteration of the inversion procedure (Fig. 15d).For these regions the variance values are higher and close to the total variance of the imposed histogram as the local conditioning with the observed ERT data is low.

[h] [h] [h] [h] [h] [h] [h] [h]
Figure 15 illustrates the integration, in a three-dimensional (3D) view, of the best-fit electrical resistivity models (Fig. 15a), as well as the pointwise mean of the models predicted in the last iteration of the inversion process (Fig. 15b).The models obtained via deterministic ERT inversion using the commercial software RES2DINV (Loke 2010) are also shown (Fig. 15c).These models were provided by the data owner and serve as a benchmark for the models predicted with the proposed ERT inversion method.The predicted models with the geostatistical inversion method show larger spatial variability, due to the stochastic sequential simulation algorithm and the imposed variogram model, and have higher coherency when interpreted together.Despite being inverted individually, there is consistency at the intersection locations; on the other hand, the results obtained with a deterministic inversion are smoother with abrupt vertical changes, which might not be geologically realistic.These abrupt variations depend on the parameterization of the inversion (e.g., vertical versus horizontal smoothing).Moreover, the integration of the deterministic solutions in a 3D view shows some resistivity spatial continuity inconsistencies, especially in the area where profiles P17 and P20 intersect.The pointwise variance model computed from the electrical resistivity models predicted during the last iteration of the geostatistical inversion is also presented in a 3D view (Fig. 15d).For the different profiles, the lowest spatial uncertainty is observed in areas where the predicted resistivity models are populated with low electrical resistivity values, while higher variability is observed in areas characterized by high electrical resistivity values.The observed ERT data for these regions tend to be smoother (i.e., with lower spatial variability) and therefore easier to match.As the observed ERT data of these regions have a higher assimilation degree during the iterative procedure the ensemble of predicted models during the last iteration has a smaller pointwise variance (i.e., spatial uncertainty).

Discussion
This study proposes an iterative geostatistical ERT inversion method based on stochastic sequential simulation and co-simulation.Information regarding the resistivity spatial continuity pattern (i.e., the variogram model) is inferred directly from the available resistivity borehole data, which might be complemented by expert knowledge.
The results obtained in both application examples show the ability of the proposed method to predict electrical resistivity models that are consistent with the recorded ERT geophysical data and with alternative deterministic inversion approaches.During the first iteration, the predicted models are already close to the true subsurface resistivity, which explains the high convergence rates.However, the iterative procedure's success and convergence rate depend on the quality of the observed ERT data, the number of boreholes, and the reliability of the estimated global variogram model.Also, in both application examples, a Wenner-Schlumberger type of acquisition array is considered.This kind of array has a good correspondence between the pseudo-resistivity section and the true spatial distribution and therefore might facilitate the convergence of the geostatistical inversion method.Tests with different acquisition geometries (e.g., dipole-dipole, multiple-gradient arrays) have shown that the performance Also, the proposed methodology is presented and illustrated with a 2D forward model; however, a 3D forward model could be used in a straightforward manner.In this case, the stochastic sequential simulation would generate sets of 3D models of subsurface electrical resistivity.
The synthetic application example is characterized by a homogeneous background versus a V-shaped contamination plume.This data set poses a challenge for geostatistical simulation methods based on two-point statistics due to the nonstationary behavior of the model parameters represented by the sharp discontinuities and the contamination plume with opposite directions.Nevertheless, in this synthetic application example, the predicted electrical resistivity models reproduce the spatial pattern of pollutant dispersion while sharing the same histogram as the reference model (Fig. 16).Nevertheless, due to the use of geostatistical simulation and co-simulation methods based on two-point statistics, the proposed inversion method struggles to predict exactly the location of the impermeable barrier and the V-shaped plume.The results agree with those obtained with the deterministic inversion.
The real case application example illustrates the method's potential with real, noisy data.The proposed method could predict spatially consistent electrical resistivity models for the different profiles from the observed ERT data.These models generated synthetic geophysical data similar to the observations while being able to reproduce the target histogram imposed for the geostatistical simulation.In this application example, the target distribution is the one retrieved with the deterministic solution, but alternative target histograms could be used (e.g., from borehole data; Fig. 15).
In the application examples shown in this paper, only the spatial uncertainty of the predicted subsurface properties is explored.However, further developments could also include uncertainty in the observed data as provided by modern ERT systems that provide an estimate of the variance of the measurements.

Conclusion
The work presented herein proposes an alternative iterative geostatistical ERT inversion method based on stochastic simulation and co-simulation.It was successfully applied to a 2D synthetic case and a set of 2D ERT profiles.This was verified by computing apparent resistivity models from the generated electrical resistivity realizations, which were locally compared against the observed one in terms of similarity coefficient.The models were constructed by selecting portions from each realization of the ensemble that showed high similarity with the observed data and then using these portions as secondary data for the next co-simulation.The ensemble of realizations generated during the last iteration of the inversion process was used to assess the uncertainty of the spatial distribution of electrical resistivity.These electrical resistivity models were characterized by variability and converged towards the areas of lower electrical resistivity values in the real case application.

Fig. 1
Fig. 1 Schematic representation of the proposed iterative geostatistical ERT inversion method

Fig. 2 a
Fig. 2 a Section of electrical resistivity reference model representing pollutant dispersion in groundwater.b Pseudo-section of reference apparent resistivity, computed considering a Wenner-Schlumberger

Fig. 3 Fig. 4
Fig. 3 Two-dimensional experimental variogram computed directly from the true electrical resistivity model and model fitted for both the a horizontal and b vertical directions.It is assumed that there is no uncertainty in the variogram model

Fig. 5
Fig.5Pseudo-sections of a reference apparent resistivity, b synthetic best-fit apparent resistivity, and c synthetic apparent resistivity given by the mean of the resistivity models after the last iteration of the inversion process.The similarity coefficient S is also shown, d for the similarity between reference and best-fit realization, and e for the similarity between the reference and the pointwise mean of all realizations.The black arrows point to small-scale differences between the reference (a) and the two synthetic estimates (b-c).w1 and w2 represent the location of the boreholes considered

Fig. 6
Fig. 6 Sections of a true electrical resistivity, b best-fit resistivity model, and c pointwise mean of the resistivity models predicted in the last iteration of the inversion process

Fig. 11
Fig. 11 Two-dimensional experimental variograms computed from the electrical resistivity data resulting from a deterministic inversion provided by the site owner: a horizontal and b vertical directions for

Fig. 13
Fig. 13 Apparent resistivity pseudo-sections of a observed and b synthetic best-fit of profile P15; c observed and d synthetic best-fit of profile P17; e observed and f synthetic best-fit of profile P19; and

Fig. 14
Fig. 14 Local similarity computed between observed and synthetic best-fit apparent resistivity for profiles a P15, b P17, c P19, and d P20

Fig. 15
Fig. 15Integration in 3D view of a best-fit resistivity model of profiles P15, P17, P19, and P20; b mean of the resistivity models predicted in the last iteration of the inversion process of profiles P15, P17, P19, and P20; c solution obtained with deterministic inversion approach of profiles P15, P17, P19 and P20; and d variance of the resistivity models generated during the last iteration of the stochastic resistivity inversion method of profiles P15, P17, P19, and P20

Table 2
Dimension of the inversion grid and global variogram model parameters used to invert profiles P15, P17, P19, and P20