Stochastic inversion of electrical resistivity changes using a Markov Chain Monte Carlo approach



[1] We describe a stochastic inversion method for mapping subsurface regions where the electrical resistivity is changing. The technique combines prior information, electrical resistance data, and forward models to produce subsurface resistivity models that are most consistent with all available data. Bayesian inference and a Metropolis simulation algorithm form the basis for this approach. Attractive features include its ability (1) to provide quantitative measures of the uncertainty of a generated estimate and (2) to allow alternative model estimates to be identified, compared, and ranked. Methods that monitor convergence and summarize important trends of the posterior distribution are introduced. Results from a physical model test and a field experiment were used to assess performance. The presented stochastic inversions provide useful estimates of the most probable location, shape, and volume of the changing region and the most likely resistivity change. The proposed method is computationally expensive, requiring the use of extensive computational resources to make its application practical.

1. Introduction

[2] A fundamental Earth sciences problem is to determine the properties of an object that we cannot directly observe. A variety of geophysical tomography techniques have been developed to provide detailed subsurface information. One such technique, electrical resistance tomography (ERT), is a relatively recent geophysical imaging scheme that provides two-dimensional (2-D) and 3-D images of resistivity that are consistent with measurements made on an array of electrodes. With the increasing availability of computer controlled multielectrode instruments and robust data inversion tools, ERT is becoming widely available. The value of ERT for monitoring dynamic subsurface processes has promoted new applications in a wide range of environments [e.g., Daily et al., 1987; Park and Van, 1991; Daily et al., 1992; Ellis and Oldenburg, 1994; Sasaki, 1994; LaBrecque et al., 1996, 1999; Binley et al., 1996; Morelli and LaBrecque, 1996; Park, 1998; Kemna et al., 2000; Slater et al., 2000].

[3] The goal of any ERT inversion method is to calculate the subsurface distribution of electrical resistivity from a large number of resistance measurements made from electrodes. A deterministic inversion procedure searches for a model (i.e., a spatially varying distribution of resistivity) that gives an acceptable fit to the data and satisfies any other prescribed constraints. A common solution minimizes an objective function consisting of a regularized, weighted least squares formulation. Typically, the search is conducted using iterative, gradient-based methods [e.g., Park and Van, 1991; Ellis and Oldenburg, 1994; LaBrecque et al., 1996].

[4] The ERT inversion problem is typically complicated by a nonlinear relationship between data and the inverted parameters, state-space dimensionality, under/over determined systems, noisy and dependent data, etc. Hence an exact inversion is rarely possible. It is common to use unrealistic simplifying assumptions to mitigate the severity of these problems when using classical optimization algorithms.

1.1. Previous Work: Stochastic Methods

[5] One alternative to the classical ERT inverse methodologies uses stochastic techniques that search for electrical resistivity models that best fit the collected electrical resistance measurements. The literature describes a variety of these methods and their application to geophysical problems. Zhang et al. [1995] suggest an inversion method that seeks to maximize a specified a posteriori probability density function of model parameters. In this case, maximizing the a posteriori density function is equivalent to minimizing the objective function in the classical inverse approach. Yang and LaBrecque [1999] proposed an alternate solution that extends Zhang et al.'s approach by allowing a more efficient estimate of the parameter covariance matrix.

[6] Mosegaard and Tarantola [1995] and Mosegaard and Sambridge [2002] describe an alternative stochastic inversion approach. The technique utilizes Markov Chain Monte Carlo (MCMC) methods, a class of importance sampling techniques that search for models that are most consistent with available data. In this approach, the inverse problem is formulated as a Bayesian inference problem. An importance sampling search algorithm is used to generate an empirical estimate of the a posteriori probability distribution based on available observations. Specifically, solutions are sampled at rates proportional to their posterior probabilities. This implies that models consistent with a priori information as well as observed data are sampled most often, while models that are incompatible with either prior information and/or observations are rarely sampled. This is the key difference between traditional Monte Carlo and MCMC: the former samples the space of possible models completely at random, while the latter moves through the models according to their posterior probabilities. This method yields an efficient sampling scheme that affords the user the flexibility to employ complex a priori information and data with non-Gaussian noise. Mosegaard and Tarantola show that this approach can be used to jointly invert disparate data types such as seismic and gravity data.

[7] Kaipio et al. [2000] describe the application of the MCMC approach to medical imaging problems using electrical impedance tomography (EIT, for the purpose of this paper, ERT and EIT are synonymous). They considered a variety of nondifferentiable priors including a minimum total variation prior, a second-order smoothness prior and an “impulse” prior that penalizes the L1 norm of the resistivity. Their approach estimates the posterior distribution of the unknown impedances conditioned on measurement data. From the posterior density, various estimates of the resistivity distribution and associated uncertainties are calculated.

[8] Andersen et al. [2001a] describe an MCMC geophysical approach for the Bayesian inversion of electrical resistivity data. They used random, polygonal models to represent the layered composition of the Earth, and demonstrate the performance of the method using field data. They analyze the posterior distribution by looking for the resistivity model that is most consistent with the data and comparing it to the estimated posterior mean model. They also estimate the variability of the transitions between Earth materials by comparing the standard deviation for each image pixel to its corresponding mean. Andersen et al. [2001b] describe another MCMC application aimed at the detection of cracks in electrically conductive media. Their approach assumes that the cracks are linear, nonintersecting and perfectly insulating. Using synthetic data, they demonstrate an updating scheme that assumes that the number of cracks is known a priori. Their approach estimates the posterior distribution of crack configurations and their associated variances.

[9] Yeh et al. [2002] describe a sequential, geostatistical ERT approach that allows inclusion of prior knowledge of general geological structures through the use of spatial covariance. Their method also uses point electrical resistivity measurements (well logs) to further constrain the solution. They use the successive linear estimator approach to find an optimal model that consists of the “conditional effective electrical conductivity” (CEEC). They define CEEC as the parameter field that agrees with electrical resistivity measurements at core sample locations and that honors electrical potential measurements. They also compute conditional variances to estimate the uncertainty associated with their optimal CEEC model.

2. Our Approach

[10] Here we describe a type of MCMC algorithm that incorporates resistance measurements, numerical forward simulators of subsurface electrical resistivity, and a priori knowledge to provide distributions of resistivity change that are likely to be present in a subsurface environment. This methodology produces resistivity models that show those subsurface plume configurations and resistivity values that are most consistent with the available data and forward models.

[11] As with all MCMC approaches, Bayesian inference is driven by an importance sampling algorithm that forms the basis for our methodology. There are two major components to the approach (refer to the flow diagram in Figure 1):

Figure 1.

Schematic diagram of the MCMC inversion process. The key differences between this approach and deterministic inversion are steps B (randomly propose inversion models) and E, F (control how the model is updated).

[12] 1. A base representation specifying the rules that the proposed resistivity models (also referred to as states) of the system must obey (step B, Figure 1). These rules are based on a priori knowledge.

[13] 2. A Markov Chain Monte Carlo (MCMC) simulation algorithm that generates samples according to the unknown posterior distribution. It uses a randomized decision rule to accept or reject the proposed states according to their consistency with the observed data (steps C–F, Figure 1).

[14] The main advantage of this approach is that it automatically identifies alternative models that are consistent with all available data and ranks them according to their posterior probabilities and associated confidences. In geophysical ERT applications, the inverse problem is substantially underconstrained and ill-posed. Thus the search for a solution that is unique and possesses a high degree of confidence is generally impossible. Hence it is wiser to consider approaches that are capable of generating alternative models and ranking them.

2.1. Base Representation: Subsurface Plumes

[15] The base representation algorithm randomly generates neighboring models utilizing rules that are based on available prior information. For example, the range of permissible resistivities is estimated using the properties of the injected fluid, and formation properties such as porosity and water saturation. We also make use of prior observations that these plumes tend to consist of one or more semicontiguous regions and relatively simple shape. Appendix A shows the specific implementation of this approach.

[16] The representation of a model space X (i.e., the individual resistivity models and their neighboring states) is critical to the overall effectiveness of the MCMC method. In those cases where generating samples from the model space is time intensive (e.g., executing forward models such as those used for three-dimensional (3-D) ERT), the model representation must be designed so that the number of degrees of freedom in the problem is well constrained. For this reason, we have chosen a categorical simulation approach where each category is associated with a discrete resistivity value. For example, if the set of categories is {A, B, C}, then AR1, BR2, CR3, where the Ri represent distinct resistivity values.

[17] A detailed outline of the algorithm that generates the resistivity models is presented in Appendix A; this algorithm corresponds to step B in Figure 1. Examples of models generated by this method are shown in Figure 2. The samples generated are called proposal states and constitute possible solutions to the inverse problem. The search for models that are consistent with the data is controlled using a Markov chain Q designed to have a “stationary” distribution equal to the prior distribution ρ(x). A stationary distribution is one in which the samples generated are representative of the true distribution because they are unaffected by the starting point of the Markov chain(s) and have explored sufficiently the distributional structure. Hence, as the Markov chain is executed, samples from ρ(x) are generated. This means that exact knowledge of the properties ρ(x) is not required because the properties are reflected by the samples from ρ(x).

Figure 2.

Examples of models generated by the base representation algorithm. Each model consists of subvolumes that can have varying size, shape, and resistivity values. The top row models show that some of the subvolumes can be separate from other subvolumes.

[18] The base representation we developed is specific to our problem domain. It only applies to resistivity models associated with liquid plumes having relatively simple shapes and contained within porous media such as those created by subsurface tank leaks or fluid injection (e.g., steam, CO2, or water floods). The algorithm assumes that the system of interest consists of a zone of changing electrical resistivity embedded within an otherwise homogeneous volume. The changing volume can be described with two or more contiguous subvolumes consisting of rectangular parallelepipeds. Each subvolume has a single resistivity value, and overlaps or is near to another subvolume. These subvolumes can have varying size, shape, and resistivity properties. This approach limits the ability of the inversion to place spurious artifacts where the data have no sensitivity. It limits the anomalies to simple shapes and significantly reduces the size of the model space.

[19] We are interested in mapping temporal changes in resistivity such as those created by a liquid penetrating the subsurface. Our approach uses a ratio of two impedance data sets in the inversion. In this method, a new data vector, dr, is formed from

equation image

where db is the data vector used as the baseline state, dt is the data vector at some time t, and fhom) is the forward solution for an arbitrarily chosen homogenous conductivity.

[20] Inversion of the new data set dr in the usual manner then results in an image that will reveal changes relative to the reference value σhom. This approach works reasonably well in situations where the contrast in background resistivity is small. Park [1998] used a similar approach except he used fheter), where σheter is the heterogeneous conductivity for the baseline case. His approach is more general and well suited for situations where the contrast in background resistivity is relatively large.

2.2. Markov Chain Monte Carlo: Theory and Methodology

[21] Our approach is a derivative of the Metropolis algorithm [Metropolis et al., 1953] as described by Mosegaard and Tarantola [1995]. It uses a Markov chain process to control the sampling of the space X of possible models. Within this framework, the solution to an inverse problem is an estimate of the posterior probability distribution defined over X. Then, for any potential solution x0X, the method will provide an estimate of the probability and confidence that state x0 is the true state of the underlying system.

[22] This MCMC approach is similar to classical inversion with the random model generator replacing the deterministic updating scheme based on a gradient search. In both cases, an initial model is chosen and responses are calculated with a forward solver. The calculated responses are compared to observed data. Finally, an updated model is chosen and the process repeats. The two approaches differ in how the updated model is chosen and the final result of the process. Specifically, MCMC produces a probability distribution defined over X, while deterministic methods produce a single or a collection of states from X that best explain the data.

[23] The inverse problem under consideration may be described as follows. Let D denote the data space, and suppose that there exists a mapping G such that

equation image

The goal is to find an x0 that corresponds to the set of observations, d0. An important point here is that X, the range of possible solutions, is limited by a priori knowledge.

[24] The sampling process can be viewed as consisting of two separate components: prior knowledge and measurements d. In section 2.1. we introduced the first component, the generation of samples that are consistent with the available a prior knowledge. We now discuss the second component: a decision process that either accepts or rejects these a priori samples according to their consistency with the measurements d (Figure 1, steps C, D, E, F). Specifically, for each visited state, forward simulators are used to predict values of measurable quantities such as electrical resistance. These predictions are generated as needed by the MCMC process, which tends to sample only the portion of solution space that is consistent with available data. These predictions are then compared to corresponding measurements to determine the likelihood L(x) that the given state xX produced the observed data. An accept/reject decision based on this likelihood is used to modify the prior sampling process. The result is a new Markov chain, R, which samples the posterior distribution, P(x). These samples provide the basis for estimating the posterior distribution and any subsequent inference concerning the true unknown state of the system.

[25] Formally, Bayes rule relates the prior and posterior distributions as

equation image

On the basis of the available prior data, the prior ρ(x) can have a distribution that is uniform, Gaussian, or unknown. We do not need to know ρ(x) explicitly because the method only requires the ability to produce samples from the distribution. The likelihood L(x) is a measure of the degree of fit between the data predicted assuming the model x and the observed data, and c is a normalizing constant. For this study, we assumed a likelihood function of the form

equation image

where N is the number of data points, d(x)pred,i is the predicted data for a given model x, d0,i is the vector of observed measurements, σi is the estimated data uncertainty, and n ≥ 1. We note that most deterministic inversions also use the term in parentheses as a measure of goodness of fit. For the results described below, we assumed that n = 2. Equation (4) assumes that the estimated data errors are uncorrelated; ERT surveys typically use the same electrodes for multiple measurements thereby increasing the probability that the data errors are correlated.

[26] The decision to accept or reject a proposed state is made on the basis of likelihood comparisons (steps E and F, Figure 1). Suppose that the current state of the Markov chain is x(T) and that a move to an adjacent state x(T+1) is proposed. If these transitions were always accepted, then the simulation would be sampling from the prior distribution ρ(x), i.e., the observed d0 would not influence the search. Instead, suppose that the decision to accept the proposed transition is made as shown by steps D, E, and F in Figure 1. Note that when the likelihood of the proposed state L(x(T+1)) is equal to or larger than that of the current state L(x(T)), the proposed transition is always accepted. If L(x(T+1)) < L(x(T)) but the two values are close to each other, the probability of acceptance is still around 1.0. For example, suppose that L(x(T)) = 10 and L(x(T+1)) = 9. As indicated by Figure 1, step E, the probability of acceptance (Paccept) will be 9/10 or 0.9. We then generate a uniformly random number RN in the range 0 to 1.0. When Paccept > RN, the transition to model x(T+1) is accepted. Note that there is a high probability of accepting x(T+1) because the odds are very high that Paccept > RN. Next, let us suppose that the model x(T+1) is much less consistent with the data such that L(x(T+1)) = 0.9. In this case, Paccept is 0.09, the odds that Paccept > RN are much smaller, and thus the odds of accepting the transition are a lot smaller. Even when L(x(T+1)) ≪ L(x(T)), Paccept is not zero. Thus this randomized rule allows a transition to a less likely state such that the process will move out of a local extremum. Theoretically, it will never get trapped in a region of locally high likelihood as long as the likelihood of the proposed state is greater than 0.0. Then, the randomized acceptance rule (step E, Figure 1) guarantees that the probability of accepting this transition will always be greater than 0.0.

[27] This is the Metropolis algorithm, the best known of the importance sampling algorithms. Metropolis et al. [1953] proved that the samples generated through this three-step process have a limiting distribution that is proportional to the desired posterior distribution P(xd), the probability of model x being the true state of nature given that d has been measured. As a result of the randomized rule in step E (Figure 1), the search tends to hover in regions of space X containing states that better fit the prior information and ERT measurements. Because of this, space X is traversed more efficiently than with traditional Monte Carlo techniques.

[28] The information contained in the ERT data determines whether the posterior distribution is a better representation of reality than the prior. When the data are “informative,” i.e., sufficiently sensitive to the characteristics of the target, the posterior distribution will be a better representation of reality. However, when the data are “uninformative” due to lack of sensitivity, measurement error, etc., then the prior and posterior distributions will be very similar, thereby indicating that the ERT data did not help to discriminate between models sampled from the prior distribution.

[29] One desirable quality of the MCMC approach is that knowledge of the posterior distribution allows the uncertainty in the generated estimate of the true unknown state to be quantified. This provides the basis for (1) the objective assessment of competing hypotheses when the available information is not sufficient to definitively identify the true system state and (2) the propagation of uncertainty in modeling results through to follow-on predictions. Sources of uncertainty such as measurement error, contradictory data, lack of sensitivity or resolution, incomplete surveys, and nonunique relationship between measurements and inverted parameters can be addressed explicitly via this approach. Moreover, problems with many secondary extrema, a nonunique inverse, and/or contradictory or sparse data are mitigated.

[30] Deterministic inversion methods are also able to address many of the problems listed above. A common approach makes use of different starting models to investigate solution uncertainty. For the models that converged, one can observe common features and gain confidence of their actual existence. This approach can be tedious to implement, and to the best of the authors' knowledge, not commonly used for ERT inversions. The MCMC approach described here performs the analysis automatically thereby sampling solution space more completely and reducing the probability of getting trapped on local extrema.

[31] There are a variety of issues that must be addressed during the implementation of the MCMC methodology. The most fundamental concern is that the Markov chain must be designed so that it has a limiting stationary distribution. For this to happen, the transition probabilities must be specified so that the process is aperiodic (state transitions are not cyclical, state sequences do not repeat) and irreducible (it is possible to move from any given state to any other).

2.3. Convergence Analysis

[32] As the Markov chain generates samples, it is important to verify that these samples are statistically representative of the posterior distribution P(xd). The chain is designed to have a long-run (i.e., stationary) distribution equal to the posterior distribution. Hence, after a sufficiently long “warm-up” period the sampling process will have forgotten its starting point and visited all of the modes (i.e., locations in space X at which a relative or absolute value occurs in the frequency distribution) of P(xd). After this point, the frequency of visits to a state will constitute a statistically reliable estimate of its posterior probability. More precisely, this means that once the chain has gone a sufficient number of steps, T0, the distribution of the generated states, x(T), at any step TT0 is unchanged and equals the posterior distribution, P(xd). We call T0 the “burn-in” period. The MCMC process begins at a particular state that is selected at random and after the burn-in period, the chain has essentially forgotten where it started. At this point, the samples image image constitute a representative random sample from P(xd) and can legitimately be used to perform posterior inference (e.g., generate estimates of the true state and their corresponding confidence).

[33] Since the burn-in period must occur before samples can be reliably collected and used, metrics have been developed to help identify T0. Several of the most widely used and effective metrics were employed during this study [Glaser, 2003].

[34] The general idea is to compare the variability within each individual chain to the variability between the separate chains by using covariance matrices. Convergence of these metrics implies that the statistical difference between the collections of samples produced by the different chains is decreasing and collectively they are beginning to exhibit the same statistical behavior. Combining this trend with the observation that their individual variations are becoming stable, provides strong evidence that all chains are sampling from the same distribution, namely, the posterior. In other words, the burn-in period has been completed and all subsequent samples are legitimately representative of the unknown posterior distribution P(xd).

[35] Gelman and Rubin [1992] describe an approach that uses multiple Markov chains to estimate the burn-in period length T0. To generate an accurate estimate, the method must address the difficulty caused by the properties and structure of P(xd) being unknown. Specifically, the posterior distribution may contain multiple modes, or likelihood peaks that the Markov chain must effectively visit in order to produce a statistically representative sample. However, since the number of significant modes is unknown a priori, we can never be certain that a single chain has explored all critical structure of P(xd). This apparent impasse is addressed through the use of multiple independent chains with individual, well-dispersed starting points. Although these chains start at different states, they share a common, but unknown, limiting distribution, P(xd). The Gelman-Rubin diagnostic effectively detects when the variability between the sample sets produced by the individual chains settles down to a value that is expected when the chains are all sampling from a common distribution (i.e., the long-run stationary distribution P(xd)). When this behavior is detected, it is likely that the burn-in process is complete.

[36] In our ERT problem, we are interested in describing contiguous subregions of specified resistivity values. For convergence analysis, each subregion is summarized by the triple Z = (z1, z2, z3), where z1 is the area, z2 is the horizontal coordinate of the centroid, and z3 is the vertical coordinate. By considering the largest contiguous subregion for each of say nine possible resistivity values, the dimensionality of the parameter of interest Z becomes p = 9 × 3 = 27. A multivariate version of the Gelman-Rubin diagnostic was used to track the behavior of several functions of the p-dimensional parameter vectors associated with the individual parallel chains for a moving and expanding window of steps (called iterations). In this study, a window is taken as a range of steps that can be characterized by a single parameter n. For example, n = 50 refers to the window of length 50 iterations ranging from iteration 51 through iteration 100, and in general, the window of size n considers each chain within the iteration sequence n + 1, n + 2, …, 2n.

[37] The quantities being tracked are correlated with the variability of Z. They include a p-dimensional matrix Wn which estimates the within chain covariances for the window n, a p-dimensional matrix Bn/n which estimates the between chain covariances for the window n and the corresponding pooled p-dimensional matrix

equation image

which estimates the covariance matrix of the posterior distribution of Z. In this last expression, m is the number of Markov chains that are running in parallel. As n increases, i.e., the window moves and expands, the influence of the starting points on the individual chains diminishes, and the following trends begin to emerge:

[38] 1. The within chain variation, summarized by the scalar det Wn, stabilizes. Typically, detWn increases as new modes are encountered by the chains in space X; detWn settles to a limiting value once all significant modes are sufficiently sampled.

[39] 2. The pooled chain variation, summarized by the scalar detV, stabilizes, a result of the combined effect of the difference between chains, characterized by B/n, becoming negligible and the within chain variation stabilizing.

[40] 3. The matrices V and W become “close” to one another. It is the “closeness” of V and W that generally indicates when the burn-in period has been achieved.

[41] The remaining challenge is how to assess the closeness of V and W. Brooks and Gelman [1998] address this issue by introducing a scalar measure of the distance between V and W:

equation image

where λ1 is the largest eigenvalue of the matrix W−1B/n. Then, as n increases, the distance between V and W diminishes, the eigenvalue λ1 decreases to 0, and Rp approaches 1. The Gelman-Rubin diagnostic, then, monitors Rp, detV, and detW, as a function of the window parameter n. For sufficiently large n, say nT0, the three conditions, Rp close to 1, detW approximately constant, and detV approximately constant, are satisfied. The nearness of Rp to 1 suggests burn-in has occurred by step T0, while stabilization of the determinants provides evidence that samples within the window starting at iteration T0 + 1 are an adequate characterization of the stationary posterior distribution.

[42] Suppose there are m chains. If the Gelman-Rubin diagnostics suggest a burn-in period of length T0, then a total of mT0 samples are discarded, and m(NT0) samples are available for analysis of properties of the posterior.

[43] In addition to the well-known Gelman-Rubin methodology, a variety of other convergence diagnostic methods have been proposed, developed and tested. These include a test of normality based on the central limit theorem [Robert et al., 1999], an examination of parameter quartiles [Raftery and Lewis, 1992a, 1992b], methods based on renewal theory [Robert, 1995], etc. In fact, the literature focused upon the assessment of MCMC convergence is quite extensive and rich. Nevertheless, at this point in time, no suite of convergence diagnostics is capable of monitoring and identifying convergence with absolute certainty.

2.4. Posterior Analysis

[44] After convergence has been verified and preburn-in models discarded, the resistivity models in the posterior distribution P(xd) can be analyzed. Our goal is to distill the relevant information in these models so that that we can infer the likely properties of the “true” resistivity model under study. The topography of P(xd) contains multiple hills whose heights are proportional to the likelihoods for each of its member resistivity models (this n-dimensional space is shown schematically in Figure 3); each point in model space represents one resistivity model. The model corresponding to the peak of each hill is commonly referred to as its mode. Multiple hills indicate that the solution to the inverse problem is nonunique, the typical case for ERT. The distribution is called multimodal when multiple hills are present and unimodal when only a single hill is present. The width of each hill indicates that there is uncertainty in the model located at the mode; this variability may be due to factors such as measurement sensitivity or measurement error.

Figure 3.

Schematic of the sector of model space included in the posterior distribution. Each grid node represents one resistivity model. Each hill represents a cluster of resistivity models having similar properties. Multiple peaks indicate that the MCMC inversion has produced nonunique results. The taller peaks identify regions containing models that are most consistent (i.e., most probable) with the observed data.

[45] This complex, multimodal structure provides a challenge when characterizing the distribution and extracting insight about the resistivity models included in P(xd). We use a clustering approach to extract this insight [Sengupta and Ramirez, 2003]. Clustering is a standard data-mining technique used to extract structure from a collection of sampled data points, in this case, sampled resistivity models. It segregates the models sampled from P(xd) into groups of models that exhibit similar properties. In our specific example, a cluster is a group of resistivity models that show similar spatial distribution of resistivity and similar resistivity values. The likelihood modes in Figure 3 represent these model clusters. The clustering process is accomplished by measuring the distance (in model space), between a model and a cluster center. A cluster's center is the model space location that best represents the central tendencies of all cluster members.

[46] When deciding whether a resistivity model should be considered a member of a particular cluster, we measure the distance (in model space) between the candidate resistivity model and the cluster center. This distance is a measure of the dissimilarity of the sample relative to the central tendencies of all the models that are already members of the cluster. A cluster's central tendencies are represented by voxel-wise distribution of resistivity values. That is, for each voxel, we calculate histograms that show how frequently each of the possible resistivity values appears in all models included in the cluster; these frequencies are normalized to lie between 0 and 1. The following example should help clarify this method.

[47] Suppose that each resistivity model contains three voxels and that the set of possible resistivities is: {10, 15, 20, 30}. Suppose further that there are 100 models in a cluster with the voxel-wise frequency distributions shown in Table 1. The frequencies for these are calculated by dividing the number of models showing a particular resistivity value by the total number of models in the cluster. Table 1 suggests that for voxel 1, there are 30 models with a resistivity of 10 (frequency is 0.3 = 30/100), 40 with a resistivity of 15, 20 with a resistivity of 20 and 10 with a resistivity of 30.

Table 1. Cluster Resistivity Frequencies (crf)

[48] We also need to calculate frequency histograms for the resistivity model in question so that they can be compared to the cluster's histograms. Suppose that the resistivity model being considered has resistivity values of (20, 15, 30). The resistivity frequencies for this model are shown in Table 2; that is, when one of the possible resistivity values is present in a voxel, the frequency for that resistivity value is set to 1.0 and frequencies for all other possible values are set to 0.0; this is repeated for all voxels.

Table 2. Model Resistivity Frequencies (mrf)

[49] We need to compare the frequencies in Tables 1 and 2 in order to calculate the model cluster “dissimilarity” MCD. For every Table 2 element where the frequency is 1.0, we subtract the corresponding element in Table 1. Thus MCD = [(1.0 − 0.2) + (1.0 − 0.2) + (1.0 − 0.1)]/3 = 0.83. Large MCD values arise when the model anomaly is located in a different part of the 3-D model, has different resistivity values, or both. When MCD approaches 0.0, the resistivity model under evaluation shows a resistivity distribution that is very similar to that of most resistivity models in the cluster.

[50] The equation for MCD can be written as

equation image

where m and n are cluster and model identification numbers, respectively, nv is the number of voxels in one resistivity model, nr is the number of possible resistivity values, mrf is the model's resistivity frequency (e.g., Table 2 values) and crf is the cluster's resistivity frequency (e.g., Table 1 values).

[51] The MCD can also be used to locate the “center state” (CS) for a given cluster. The CS is that resistivity model showing the minimum MCD; it is also the model that is closest to the cluster's “center of mass” and should be the one that best represents the cluster members. The CS and the mode refer to the same model when the cluster members are distributed symmetrically about the mode.

[52] The algorithm we used to perform the clustering analysis is called the dynamic k means algorithm described in detail in Appendix B. The flow diagram for this procedure is presented in Figure 4.

Figure 4.

Flow diagram for the “dynamic k means” clustering algorithm used for this work.

[53] Once partitioning of the models is complete, the voxel-wise average resistivity, VARi, is used to provide a reasonable estimate of the trends exhibited by most models in that cluster. Consider a column vector u with nv components, the components containing the voxel resistivities in a 3-D resistivity model. If there are N column vectors corresponding to all the resistivity models in a given cluster, then

equation image

This is repeated for all components (i = 1, 2, …, nv). For non-ERT applications, one may choose to use an arithmetic mean instead of the geometric mean indicated by equation (8). The geometric mean is preferred because of the wide range in electrical resistivities.

3. Results and Discussion

3.1. Physical Model Results

[54] We have used physical models where we know the exact properties of the target in order to evaluate the performance of the MCMC approach. The physical model consisted of various targets immersed in a tank filled with water. The model included 4 vertical electrode arrays, each having 15 electrodes (refer to Figure 5). The arrays were submerged in water and a variety of solid and porous (sand-lead mixture encased in a nylon mesh) targets having different electrical resistivities were inserted at various locations between the electrodes. The water resistivity was 16 ohm m. We will present results from the sand-lead mixture that had a resistivity of about 40 ohm m.

Figure 5.

Schematic representation of the physical model setup. Four vertical electrode arrays were immersed in a fiberglass water tank. Various objects were inserted between the arrays at a variety of locations and ERT data collected.

[55] Issues such as the accuracy of the inverted location, shape and change magnitude of the inversions were evaluated. Uncertainty arises from the inherent errors (measurement and modeling) and the nonunique relationship between inverted parameters and measurements.

[56] We first needed to determine if the inversion had converged by using the convergence diagnostics described earlier. The Gelman-Rubin diagnostic is a metric indicating when T0 (burn-in period length, stable value near 1.0) has been reached asymptotically. For the physical model results discussed below, the burn-in period ended around iteration T0 = 700. Also, stable values near 1.0 suggested that samples after T0 had succeeded in visiting all the modes in P(xdo). This diagnostic suggested that the Markov chains converged to the limiting posterior distribution P(xdo), and that we could reliably employ these posterior samples to estimate the parameter(s) of interest.

[57] Figure 6 shows the stochastic inversion results for the sand-lead physical model. Figure 6 (left) is a vertical section of the actual target showing its shape and location. Each image consists of 16 × 16 × 33 voxels representing a total volume of 6.8 × 10−2 m3; the volume of the target is 2.5 × 10−5 m3. The model has a relatively low-contrast (resistivity only 2.5 times as high as the surrounding water). For these inversions, we assumed that the target could be sufficiently described by two contiguous or overlapping parallelepipeds and that the set of possible resistivity ratio values was {0.5, 1.0, 2.0, 4.0, 8.0, 16.0}. A value of 1.0 indicates that the target's resistivity equals that of the surrounding water, and values above 1.0 indicate that the target's resistivity is larger than the surrounding water. The inverse process searched for the most likely location, size, shape and contrast of the changing region.

Figure 6.

Clustering analysis results for the case of the sand-lead target. The 3-D block shown corresponds to the volume enclosed by the dashed lines in Figure 5. (top left) Vertical slice through the target. The vertical slices were placed where the maximum resistivity ratio is observed. (top right) Voxel-wise average resistivity obtained when a unimodal posterior distribution is assumed. (bottom) Average resistivity ratio obtained when a multimodal distribution is assumed. The three most probable clusters and their corresponding frequencies are shown.

[58] To calculate the inversion results shown, we used all postburn-in models (all x(T) such that T > 700, a total of about 4300 models). Figure 6 (top right) shows the VAR for all the posterior models assuming that all models came from a unimodal distribution (single cluster); i.e., no attempt was made to segregate the resistivity models in the posterior distribution into groups of models having similar spatial distribution of resistivity and similar resistivity values. The vertical slice shown through each model was positioned to show the maximum resistivity value present. The voxel-wise average results show similar shape and location as the actual target. The resistivity ratio magnitude lies in the range of 2.0–2.9, while the actual target value is 2.5.

[59] We next considered the possibility that posterior models came from a multimodal distribution. Using the dynamic K means technique described earlier, u = 10, and td = 0.04, we obtained the results shown in Figure 6 (bottom) (where td is the user-defined, model space distance; when the minimum of all distances between a model and existing cluster centers exceeds td, a new cluster is created; additional details are presented in Appendix B). The analysis produced 10 clusters of models; the three most probable ones (showing the highest frequency) are displayed in Figure 6 (bottom). The results shown are cluster centroids (VAR of all the models in a given cluster). The frequency values have been normalized relative to the number of posterior distribution samples (4302) to obtain relative frequency values for each cluster. Figure 6 (bottom left) is the centroid for the most frequent cluster: four out of every five samples come from this part of the state space. Note that the location and size of the cluster centroid are similar to the target. The maximum value resistivity ratio magnitude is about 2.48 while the target value is 2.51. Figures 6 (middle) and 6 (right bottom) represent centroids for low-frequency clusters. Figure 6 (middle) shows two closely spaced anomalies that straddle the target's elevation and a resistivity ratio value (16.0) that is substantially larger than the target. Figure 6 (bottom right) also shows an anomaly shape that is significantly different from the target and has a maximum resistivity ratio value (16.0) that is substantially larger than the target.

[60] We propose that the results shown in Figure 6 illustrate the value of clustering analysis. When the cluster analysis accounts for several likelihood hills in the posterior, the resistivity ratio for the most frequent cluster is about 2.48, close to the target value of 2.51. When a unimodal distribution is assumed, a poorer match to the target is observed because all samples have lumped into a single cluster. The clustering analysis segregates posterior samples with lower and higher values into separate clusters, thereby improving the accuracy of the resistivity value corresponding to the most frequent cluster. Similar comments apply to the size, location and shape of the cluster anomalies.

3.2. Field Results: Tank Leak Detection

[61] Electrical resistance data were collected during a field experiment that simulated leakage from a large metallic tank. For testing, an electrical tracer (saline solution) was used instead of the real contaminant to preserve the environmental quality of the test site. The test site used for this work is part of the 200 east area in the Hanford site, located near Richland, Washington (additional details about the test site and testing are provided by Barnett et al. [2002, 2003]). The near-surface sediments at the test site consist primarily of fine to coarse-grained sand with an average fractional porosity of 0.25 to 0.30. The field experiments were performed under a 15.2 m diameter steel tank mockup. Figure 7 shows the layout at the experimental site. This empty steel tank contained several built-in spill points (only the one used is shown). The bottom of the tank is located 1.5 m below ground surface. Sixteen boreholes with eight electrodes in each surrounded the tank. The electrodes were spaced every 1.52 m between the surface and 10.7 m depth. The water table is located tens of meters below the bottom electrode. The diametrical distance between boreholes was 20.7 m.

Figure 7.

Schematic layout of the leaking tank site. Hypersaline brine solution was released from a point near the center of the tank's bottom. Sixteen vertical arrays of electrodes were used to monitor the infiltration process.

[62] Hypersaline solution was released from a point near the tank's center over a 52 hour period. The liquid consisted of a sodium-thiosulfate solution (36 wt%) with a conductivity of about 5 S/m (about 1.5 times that of seawater). This hypersaline fluid has similar electrical conductivity and density to the real Hanford tank liquids. ERT surveys were made before, during and after the brine release using a dipole-dipole approach.

[63] Tank leak results are shown in Figure 8 for the case where 2160 L had been released. Dimensions of the 3-D image block are 28 m (along each horizontal axis), and 13 m (height). Each model consists of 22 × 22 × 17 voxels. For these inversions we assumed that the target could be sufficiently described by six contiguous or overlapping parallelepipeds. Their shape, location and resistivity contrast were allowed to vary. The set of possible resistivity ratio values was {1.0, 0.95, 0.90, …, 0.01}, for a total of 20 possible values. The posterior samples were analyzed using the “dynamic k means” algorithm described earlier with td = 0.04 and u = 20. The clustering was done using all postburn-in models (all x(T) with T > T0, a total of 3802 models; T0 ≈ 810). Figure 8 (left) displays the three most likely cluster centroids with frequencies ranging from 0.12 to 0.38. These three clusters encompass 68% of all posterior distribution samples. All of them suggest a roughly vertical anomaly directly below the release point to depths ranging from 8 to 10 m.

Figure 8.

Clustering analysis results corresponding to the tank brine release experiment. The 3-D block shown here corresponds to the 3-D block located beneath the tank in Figure 7. The 3-D block shown here corresponds to the 3-D block under the tank shown in Figure 5. (left) Voxel-wise average resistivity ratio for the top three most probable clusters. (right) Center state for the three most probable clusters.

[64] Figure 8 (left) shows the voxel-wise mean resistivity ratio for the three most probable clusters. The frequency values shown have been normalized relative to the number of posterior distribution samples. Clusters A and B suggest that a strong vertical anomaly exists just east of the release point. Cluster A suggests the possibility of liquid ponding to the west of the release point at a depth of about 8m. Likewise, cluster C suggests a wider but somewhat weaker anomaly. An analyst could consider cluster A as the most likely thereby inferring a pillar-shaped vertical invasion zone with possible liquid ponding to the west of the release point at 8m depth. Considering clusters A and B, the analyst may also consider the possibility that the bulk of the contamination is located below and to the east of the release point. This approach provides the analyst alternative models with model A being the most likely, and model C the least likely.

[65] This analysis can be used to evaluate the propagation of uncertainty due to measurement error or due to lack of sensitivity or resolution. For example, the mean resistivity ratios (Figure 8, left) for all three clusters show significant variability directly below the release point. The highly conducting metal tank walls create a region of diminished sensitivity below the tank that is most severe just below the release point. We suggest that this diminished sensitivity is responsible for the variability observed. However, the method cannot distinguish between variability due to data misfit and variability due to lack of sensitivity.

[66] Figure 8 also shows the center states for the clusters shown. Within a given cluster, the differences between the voxel-wise means (Figure 8, left) and the center states (Figure 8, right), provide the analyst some measure of variability (uncertainty) within the cluster. These center states are those posterior samples having the smallest MCD. We can use the center state as a way to evaluate a cluster's central tendencies. These center states are unaffected by the variability of such properties within a cluster; clearly, this variability affects the mean values shown in Figure 8 (left). The center states for clusters A and B show similar location and shape as the cluster mean image. However, the mean for cluster A suggests the possibility of liquid ponding to the west of the release point at a depth of about 8 m, whereas the cluster center state does not show such features. This means that the possible ponding of liquid is not part of the central tendency for cluster A, but is a part of some of the states around the cluster center. The center state for cluster B suggests a stronger resistivity change (0.03) than the values near 0.2 indicated by the cluster's mean. This approach allows an analyst to consider alternative models that are consistent with the data to varying degrees. The frequency information determines how well each of the alternative models represents the available data.

[67] A drill-back program that would have independently mapped the plume's characteristics was not carried out because of lack of funds. Instead, we are forced to use circumstantial evidence to evaluate the results. First, we visually compare the stochastic inversion to a deterministic tomograph calculated by the inversion algorithm described by LaBrecque et al. [1999]. We have confidence that this algorithm produces reasonable results because it has been used successfully in field applications [e.g., LaBrecque et al., 1999; Morelli and LaBrecque, 1996; Daily and Ramirez, 2000] as well as in controlled physical model experiments where the target characteristics are known exactly [Ramirez et al., 2003]. Figure 9 shows the comparison of the two results. The stochastic result (Figure 9, left) shows the voxel-wise mean for the most likely cluster (cluster A, Figure 8, top left). Figure 9 (bottom) shows isosurface views, where the white bar over the color bar indicates the range of transparent values used to render the isosurfaces.

Figure 9.

Comparison of (left) stochastic and (right) deterministic inversions. The 3-D block shown here corresponds to the 3-D block under the tank in Figure 5. (top) Series of vertical slices oriented parallel to north-south line. (bottom) Isosurfaces. For the isosurfaces, all values >0.9 are transparent.

[68] Figure 9 shows similarities between the classical inversion and MCMC results thereby suggesting that the stochastic results are reasonable. The location of the changing region is approximately the same. Also, both methods show that the zone of change extends more to the east of the release point than to the west. There are also differences between the images. The stochastic image shows larger resistivity changes reaching values of 0.1 (i.e., the invaded soil has 1/10 the resistivity of the prerelease value) and a smaller volume. The stochastic image also suggests that brine may be pooling on the western side of the image. In the absence of independent ground truth data we cannot establish which of the two images is closer to reality.

[69] We can compare the two results on the basis of how well either fits the resistivity measurements. We calculated the forward solution for the deterministic and stochastic inversions shown in Figure 9. Then, we calculated the root mean squared differences (RMS) between the forward solutions and the measurements. The RMS for the deterministic result is 1.6 × 10−7, substantially better than the MCMC result (1.7 × 10−5). We speculate that the deterministic result fits the data better because the deterministic inversion allows a much wider range of resistivity values and anomaly shapes. We can also compare the volume and average resistivities of the inversion anomalies to invaded zone estimates based on Archie's equation and independent observations of released brine resistivity (0.11 ohm m) and volume (2160 L). We assume that average porosity is about 0.25, the resistivity of the prerelease pore water is about 5 ohm m, and that the brine release increases the saturation from 0.4 to 0.8. We also assume that the brine displaces part of the original pore water so that pore water resistivity after release has decreased to 2.5 ohm m. Given these assumptions, brine would change the bulk resistivity from 316 ohm m to about 40 ohm m (resistivity ratio of 0.13) within a 10.8 m3 volume of soil. For the deterministic inversion anomaly, the geometric mean resistivity ratio is 0.96 within an anomaly volume of 5511 m3. For the MCMC inversion, the geometric mean resistivity ratio is 0.21 within an anomaly volume of 1387 m3. If instead of the cluster mean, we choose the center state for cluster 4 (Figure 8, top right), then the resistivity ratio is 0.19 within an anomaly volume of 391 m3, the closest to the estimate.

[70] Clearly, both methods produce anomaly volumes that are grossly exaggerated probably because of the intrinsically low resolution of the ERT measurements and to decreased measurement sensitivity caused by shunting of electrical current through the tank's metal walls. It is also possible that fluid flow maybe occurring via capillary action along grain surfaces, thereby affecting a much larger volume. The comparison suggests that the MCMC inversion produced an anomaly that is more compact and exhibits a larger resistivity change than the deterministic result. These differences are at least partly due to the regularization used by the deterministic algorithm that penalizes models with large “roughness” (the inverse of smoothness) [Park and Van, 1991; Shima, 1992; Ellis and Oldenburg, 1994; Sasaki, 1994; LaBrecque et al., 1996; Morelli and LaBrecque, 1996]. It is well known that this approach tends to produce models that have reduced contrast and exaggerated extent.

[71] In this work we have compared an MCMC method that is constrained by size, location, and resistivity limits with a deterministic method that is constrained by model roughness. The ideal comparison between the deterministic and stochastic methods would compare the results using the same set of constraints. We suggest that such a comparison be conducted as part of future research.

[72] Figure 10 shows further evidence that the stochastic results are reasonable. Figure 10 shows a sequence of MCMC inversions as the volume of released fluid increases from 340 to 2160 L. Along a given column, released volume increases from Figure 10 (top) to Figure 10 (bottom). Figure 9 (right) shows an isosurface view where resistivity ratios from 0.9 to 1.0 are transparent. As expected, the voxel-wise average ratio decreases (i.e., the soil becomes more electrically conducting) from about 0.95 (Figure 10 (top) released volume 0.34 m3) to about 0.2 and below (Figure 10 (bottom) released volume of 2.16 m3). Also, the vertical and horizontal extent of the anomaly grows with increasing released volume. The sequence shown suggests that the plume represented by these inversions is behaving as expected: anomaly volume grows and the resistivity ratio decreases as released brine volume increases.

Figure 10.

Series of stochastic inversions corresponding to the brine release sequence; released brine volume increases from top to bottom. (left) Voxel-wise average resistivity ratio and (right) same results as isosurfaces. For the isosurfaces, all values >0.9 are transparent.

[73] Last, we suggest that the stochastic results in Figures 810 are consistent with plume characteristics determined independently by Ward and Gee [2001] and Gee and Ward [2001]. Their field tests were conducted at a site located a few hundred yards from the Mock Tank leak test facility where our results were obtained; both sites are known to have similar geology. Their field tests, laboratory tests, and numerical simulations of the movement of hypersaline solutions through the vadose zone suggest that these plumes move along finger-like, vertical preferential flow paths because of the much higher density of the hypersaline solution. Their results also show that the plumes tend to be more compact, move deeper and show less lateral spreading than low ionic strength plumes. The anomalies shown in Figures 810 show similar characteristics to those observed by Ward and Gee and Gee and Ward.

[74] In summary, the similarity between the stochastic and deterministic results, the increasing stochastic anomaly volume as the released brine volume increases, and the consistency between stochastic anomaly characteristics and independent a priori expectations of the plume characteristics, leads us to believe that the MCMC results provide reasonable representations of reality.

3.3. Computational Expense: Parallel Computing

[75] The MCMC method we have used is computationally expensive. For example, a tank leak problem involving 28,800 voxels and 128 transmitting electrodes and a typical work station with one CPU, will require about 45 days to accumulate a sufficient number of posterior samples (about 4000). Almost all of that time is used to solve the forward problem. In comparison, a deterministic inversion of the same problem using the same workstation may take 8–12 hours.

[76] Clearly, parallel computation of the forward problem is required to make the MCMC approach practical. We have parallelized the problem in two ways: (1) Individual Markov chains are run on separate processors. (2) The computational load for each chain is further distributed among multiple processors by computing the potential field due to each transmitter electrode on separate processors. The potential fields are then postprocessed to compute the transfer resistances associated with individual readings. When 128 processors are used, this approach reduces the processing time to about 12 hours. Processing times of 12 hours are acceptable for many real-life applications of the MCMC approach. Clearly, approach 2 could also be used to accelerate a deterministic run.

[77] Step size (i.e., distance between neighboring models in space X) controls the characteristic changes one allows in model space. Step size has a large impact on computational expense; if the resolution is too fine or involves a high dimensional state vector, the convergence may be slowed beyond practical limits. If the step size is too small, movement through the state space will be slow, and it will take longer to move past local extrema. When the step size is too large, increased rejection ratios (number of states rejected/total number of states evaluated) are likely thereby slowing the convergence rate. By trial and error, we discovered that a randomized step size that sometimes took smaller or larger steps provided a reasonable solution to this dilemma.

[78] Closely related to the step size is the topography of the likelihood surface under analysis. When the surface is exceptionally steep, the process will be slowed because of high rejection rates as the process attempts to move off a steep peak in likelihood space and ends up proposing states having much smaller relative likelihoods that are almost always rejected. This slow mixing process (i.e., how efficiently the process moves through space X) can be mitigated in several ways but that discussion is beyond the scope of this manuscript. Finally, the choice of the prior distribution may also significantly impact convergence because it controls the proposal of candidate states. In general, the closer the prior is to the posterior, the faster the process converges to P(xdo).

4. Summary and Conclusions

[79] We have discussed a stochastic methodology for the inversion of changing subsurface electrical resistivity data. This method is based on Bayesian inference and is implemented via an MCMC algorithm. The inversion of electrical resistivity data is an ill-posed problem requiring regularization. Our approach makes use of prior information to sufficiently reduce the size of the space of feasible solutions in order to mitigate ill-posedness. The resistivity models consist of multiple blocky regions of resistivity change embedded within an unchanging volume. Additional information can include the sense of the change (increasing or decreasing resistivity), upper/lower bounds for the volume of the changing region, resistivity change magnitude and spatial relationships of the regions (e.g., requiring partially overlapping or contiguous blocks).

[80] A key strength of MCMC is that solutions are sampled at a rate proportional to their consistency with available data. Hence models that are most consistent with available data observations are sampled most often, while models that are incompatible with either prior information and/or observations are rarely sampled. As a result, the frequency of models in the posterior distribution can be used to determine the probability that a given model is the best explanation for the available data. The approach can be used to identify competing models when the available information is not sufficient to definitively identify a single optimal model. Another strength is that it can be used to jointly invert disparate data types such as seismic and gravity as shown by Mosegard and Tarantola [1995].

[81] We view the MCMC and the deterministic inversion methodologies as complementary approaches. The MCMC approach is similar to classical inversion with the random model generator replacing the deterministic updating scheme based on a gradient search. In both cases, an initial model is chosen and responses are calculated with forward solver. The calculated responses are compared to observed data. Finally, an updated model is chosen and the process repeats. The two approaches differ in how the updated model is chosen.

[82] The deterministic method is likely to be the preferred method when fast inversion times are required, and when the regularization scheme produces sufficiently accurate models. The MCMC method may be most useful when inverting problems with many secondary extrema, when explicit estimates of solution uncertainty are required, and when alternative models are desired, ranked according to their consistency with available data.

Appendix A:: Base Representation Algorithm: Generate Proposal Models

A1. Introduction

[83] Here we describe the algorithm used to generate the resistivity models. The algorithm assumes that the resistivity model consists of a mass of changing electrical resistivity embedded within an otherwise homogeneous volume. The changing volume is composed of two or more subvolumes consisting of rectangular parallelepipeds. Each subvolume has a single resistivity value selected from a user-specified set of values; other subvolumes can have different resistivity values. A subvolume has to overlap, be contiguous to or be near another subvolume. These subvolumes can have varying size, shape, and resistivity properties.

[84] Let xT and xT+1 be two adjacent resistivity models in space X; xT is the current model in the Markov chain and xT+1 will be the new model. Let nvmax be the maximum number of subvolumes that is allowed.

A2. Algorithm

[85] For each subvolume making up the mass of resistivity, propose new size, new location, and new resistivity values for xT+1 by perturbing the xT values.

[86] Step 1 is for each subvolume sv (1 ≤ sv ≤ nvmax) in xT+1:

[87] Step 1a proposes the size in elements along the X, Y, Z directions by randomly perturbing the size of model xT.

[88] If proposed size is outside acceptable range, flag the proposed size.

[89] Step 1b proposes a new location for the each subvolume by randomly perturbing the location for xT; do this for the X, Y, Z directions.

[90] Flag the proposed change if part or all of the subvolume is outside the allowable region.

[91] Step 1c is to randomly move through the set of allowable resistivity values to select a new resistivity values for each subvolume. Suppose that the set of possible resistivity values of {1.0, 2.0, 8.0, 16.0, 32.0} and let the resistivity value for xT = 8.0 (position 3 on the set).

[92] Choose a random integer in the range −2 to 2 (or some other user-specified range) and add it to the position number corresponding to xT. The resistivity value for xT+1 will be the value associated with this new position. Thus, if the random integer = 1, the resistivity value for xT+1 is the value associated with position 4, 16.0.

[93] Check that proposed set element is one of the positions within the set; flag it if it is not.

[94] Step 1d checks whether the proposed new size, location or resistivity value are outside the range of permissible values; if so, go back to step 1a and try again.

[95] Step 2 checks that each subvolume is contiguous, overlaps or is acceptably close to at least one other subvolume, do this check along the X, Y, Z directions.

[96] If each subvolume is not contiguous (or acceptably close) to at least one other sub volume along all three principal directions go back to step 1.

[97] The save_previous_function step saves the values of the xT model (size, location, resistivity model). This function is called before proposing the values for the xT+1 model. Do this for all subvolumes.

[98] The reset_function is used to reset the Markov chain back to the xT model whenever the proposed xT+1 model is rejected by the Metropolis algorithm. This function sets all the values for proposed model xT+1 equal to those of xT by using the values saved by “save_previous_function”.

Appendix B:: Dynamic K Means Clustering

[99] We describe the clustering algorithm used (refer to Figure 4). A variation of the K means [Jain and Dubes, 1988] algorithm is obtained when starting from an initial value of K = 1 clusters, we allow K to grow dynamically as new resistivity models are considered for assignment to an existing cluster [Pao, 1989]. A new cluster is formed when a sample due for assignment to a cluster is considered as being too “far” (i.e., has too large an MCD as defined by equation (7)) from any existing cluster. Clearly, a new parameter td, the threshold distance has to be introduced at this stage. Then, when the minimum of all distances between the model and existing cluster centers exceeds td, we create a new cluster with the sample under consideration as its first member. The process is repeated iteratively until a specified stability criterion is satisfied.

[100] Suppose that K is the number of clusters, u is the maximum number of clusters allowed, and td is the threshold distance.

[101] In stage 1 (cluster growing), initially, set K = 1 with cluster 1 containing the first model and the cluster center located at the first model. At any given stage, repeat steps 1 and 2 where,

[102] 1. Get a new model and compute its MCDs from the existing cluster centers and their minimum(MCD) and locate a cluster C say, where this minimum occurs.

[103] 2. If minimum(MCD) is greater than td, then create a new cluster with the new model as the only cluster member and itself as the new cluster center.

[104] Otherwise, assign the new model to cluster C and update the cluster center. Continue this process until either the maximum number of clusters has been attained or all of the models are assigned to one of the generated clusters.

[105] In stage 2 (reclustering), at the end of the first stage, let there be Ku clusters generated.

[106] Repeat steps 3 and 4 until a stability criterion S (defined below) is satisfied.

[107] 3. Replace the existing partition with a new one by assigning a group to each member based on the nearest cluster principle (minimum MCD).

[108] 4. Recompute the cluster centers by determining the center state for the new collection of members in each cluster.

[109] Tests of this algorithm show that stage 2 insures that the clustering results are independent of the order in which the models are introduced. The stability criterion S can be introduced in various ways. In our algorithm, S is reached when the updated cluster centers all remain within a small preassigned distance from the respective cluster centers computed in the previous step.

[110] Two clusters may be merged to form one whenever the maximum of the intracluster distances within each cluster is smaller than the intercluster distance between the two clusters. Suppose that the “diameter” of a cluster is defined as the maximum of the intracluster distances between models within a cluster. Let us also define the “intercluster distance” between two clusters as follows. Consider all pairs of models that can be formed by choosing one model from each of the two clusters. The largest distance between models within all such pairs is then defined as the intercluster distance. In order that two clusters may be merged to form into one, we require that the sum of their diameters be less than the intercluster distance between the two. Moreover, this final merging step does not require any distributional assumptions or the specification of additional parameters to be implemented. The new cluster center and the corresponding cluster frequency are computed and cluster memberships for the merged clusters reassigned.


[111] This work was funded by the Laboratory Directed Research and Development Program at Lawrence Livermore National Laboratory. Robin Newmark (LLNL) assisted in getting funding for the project and getting it started. We value the contributions of JGR Associate Editor Stephen Park and two anonymous reviewers; their insight and recommendations made the final version of the manuscript much better than the original. We also thank Andrew Binley, Lancaster University (UK), who served as an informal reviewer. DOE's Office of River Protection, CH2M Hill Hanford Group sponsored the leaky tank experiment. The experiment was designed and managed by G. Gee, B. Barnett, and M. Sweeney, Pacific Northwest National Laboratory. The physical model data were collected by Stan Martins (LLNL) as part of a Laboratory Directed Research and Development project led by Charles Carrigan (LLNL). This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under contract W-7405-ENG-48.