Journal of Applied Geophysics 96 (2013) 38–54

Contents lists available at SciVerse ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

A Bayesian trans-dimensional approach for the fusion of multiple geophysical datasets Arash JafarGandomi a,b, Andrew Binley a,⁎ a b

Lancaster Environment Centre, Lancaster University, Lancaster LA1 4YQ, UK School of GeoSciences, University of Edinburgh, Edinburgh EH9 3JW, UK

a r t i c l e

i n f o

Article history: Received 16 October 2012 Accepted 3 June 2013 Available online 13 June 2013 Keywords: Inverse theory Data fusion Markov-chain Monte-Carlo Trans-dimensional Electrical resistivity tomography

a b s t r a c t We propose a Bayesian fusion approach to integrate multiple geophysical datasets with different coverage and sensitivity. The fusion strategy is based on the capability of various geophysical methods to provide enough resolution to identify either subsurface material parameters or subsurface structure, or both. We focus on electrical resistivity as the target material parameter and electrical resistivity tomography (ERT), electromagnetic induction (EMI), and ground penetrating radar (GPR) as the set of geophysical methods. However, extending the approach to different sets of geophysical parameters and methods is straightforward. Different geophysical datasets are entered into a trans-dimensional Markov chain Monte Carlo (McMC) search-based joint inversion algorithm. The trans-dimensional property of the McMC algorithm allows dynamic parameterisation of the model space, which in turn helps to avoid bias of the postinversion results towards a particular model. Given that we are attempting to develop an approach that has practical potential, we discretize the subsurface into an array of one-dimensional earth-models. Accordingly, the ERT data that are collected by using two-dimensional acquisition geometry are re-casted to a set of equivalent vertical electric soundings. Different data are inverted either individually or jointly to estimate one-dimensional subsurface models at discrete locations. We use Shannon's information measure to quantify the information obtained from the inversion of different combinations of geophysical datasets. Information from multiple methods is brought together via introducing joint likelihood function and/or constraining the prior information. A Bayesian maximum entropy approach is used for spatial fusion of spatially dispersed estimated one-dimensional models and mapping of the target parameter. We illustrate the approach with a synthetic dataset and then apply it to a field dataset. We show that the proposed fusion strategy is successful not only in enhancing the subsurface information but also as a survey design tool to identify the appropriate combination of the geophysical tools and show whether application of an individual method for further investigation of a specific site is beneficial. © 2013 Elsevier B.V. All rights reserved.

1. Introduction A variety of near surface geophysical methods have been recognised as potentially useful tools for efficient and reliable site characterization. However, each method is constrained by its achievable spatial resolution and coverage, ability to delineate the inherent heterogeneity of the subsurface and practical implementation over a site. Consequently, there is growing interest in the assimilation of data from multiple geophysical methods. The applicability of different geophysical methods depends on: (1) the nature of the subsurface and the target of interest (e.g. contamination, hydrology, etc.); (2) the geography of the site, since this limits the practical applicability of geophysical methods; and (3) sources of noise that may exist at the site.

⁎ Corresponding author. Tel.: +44 1524 593927. E-mail addresses: [email protected] (A. JafarGandomi), [email protected] (A. Binley). 0926-9851/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jappgeo.2013.06.004

Fusion of geophysical datasets obtained from different geophysical methods has been an active subject of recent research. Multiple geophysical methods have been successfully used to characterize hydrologic systems (Doetsch et al., 2010a,b; Linde and Doetsch, 2010; Looms et al., 2008; McClymont et al., 2011; Titov et al., 2005); to reduce the subsurface uncertainties associated with hydrocarbon and CO2 reservoirs (Chen et al., 2007; Commer and Newman, 2009; Hu et al., 2009; JafarGandomi and Curtis, 2012; Rovetta et al., 2008); to improve investigations on archaeological sites (Kvamme, 2006, 2007), and many other applications. One data fusion approach is to invert different datasets individually to obtain separate subsurface images and then combine these images and interpret them simultaneously (e.g., McClymont et al., 2011). However, most attempts to fuse multiple datasets in a formal procedure have been based on joint inversion, either within a probabilistic or deterministic framework, in an area of investigation that has been surveyed by more than one geophysical tool. The structural approach formulated by Gallardo and Meju (2004), offers great potential for such

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

application, although petrophysically-based techniques have also been successfully applied (e.g., Hoversten et al., 2006). While these methods can provide improved characterisation, they have been limited to investigations where each technique is applied to exactly the same subsurface region. In many practical applications of near surface geophysical methods, however, this is a seriously limiting constraint. Often, geophysical methods have different spatial coverage at a site due to constraining factors in each method. Whilst, traditionally the results from these multiple approaches are fused in overall interpretation using expert knowledge, there is a need for a more formal, objective approach for such data fusion. Furthermore, as we seek to improve our knowledge of a site, it may be instructive to provide a means of evaluation of the information content in each potentially applicable technique. The main objective here is thus to develop a consistent fusion approach for multiple geophysical datasets and to assess the applicability of different geophysical methods, individually and in combination with each other, and to reduce uncertainties in geophysical parameter estimates. Reducing uncertainty requires elaborate data gathering, processing and inversion workflows. Optimal experimental and data processing designs may be used to improve these uncertainties (e.g., Curtis, 2004a,b; Maurer et al., 2010). The choice of method used to design experiments depends greatly on how one can measure the information about model parameters that is expected to be gained from the data. Here we draw from experimental design and information theory to develop methods to estimate information about the variation of electrical resistivity, using this as a key geophysical parameter that assist in the characterization of a site. This, in turn, identifies key contributions to information and allows an appropriate selection of geophysical technique(s) to reduce overall uncertainty of the subsurface. Over the past two decades a number of algorithms have been developed for two-dimensional (2D) inversion of electrical resistivity data (e.g., Binley and Kemna, 2005; Oldenburg and Li, 1994). However, in two-dimensional inversion the layer boundaries are usually smeared and, consequently, the earth models produced are smoothed. Furthermore, incorporation of prior information and joint inversion of multiple datasets retrieved with different methods is challenging. There have been efforts to combine one-dimensional (1D) and 2D calculations to increase the computation efficiency (Smith and Booker, 1991) and to improve spatial resolution (Gyulai and Ormos, 1999; Santos, 2004). Auken et al. (2005) proposed a piecewise 1D laterally constrained inversion (LCI) to produce laterally smooth models with sharp layer interfaces. Brodie and Sambridge (2006) developed a quasi-three-dimensional (3D) layered inversion methodology, which is specially designed for airborne electromagnetic data. They use bicubic B-splines to invert a 3D grid (using a 1D forward solution) for a combination of layered-earth parameters. Viezzoli et al. (2008) expand the concept of LCI to spatially constrained inversion (SCI), which operates on 3D datasets and is similar to the quasi-3D inversion methodology of Brodie and Sambridge (2006). Although the above algorithms are efficient, they do not address the appraisal of the estimated models in the inverse problem (i.e., assessing representativeness of the estimated models with respect to the true model). In this work, we apply a trans-dimensional Markov chain Monte Carlo (McMC) approach to address non-linearity, non-uniqueness and uncertainties in the inversion of multiple geophysical datasets. Our fusion strategy is based on the capability of multiple geophysical methods to provide enough resolution to identify subsurface material parameters and structure. Recently, Minsley (2011) outlined a trans-dimensional McMC algorithm for spatial mapping from airborne or ground-based electromagnetic conductivity surveys. We adopt a similar philosophy here but consider data from multiple geophysical modalities, in addition to examining the information content of the post-inversion result (posterior

39

probability distribution). We develop quasi-2D and 3D approaches, which consist of re-casting the subsurface model to a set of locally 1D earth-models. To estimate parameters of each 1D model, the 2D and 3D data are either re-casted to a set of equivalent 1D data (e.g., trans-dimensional 2D DC resistivity data to a set of vertical electric soundings (VESs)) or interpreted separately to obtain corresponding 1D earth-models (e.g., ground penetrating radar (GPR) data), which provide prior information for the joint inversion. Selection of the action to be taken depends on the data type. These estimated 1D earth-models are then stitched together to construct 2D images or 3D blocks of the subsurface model. The uneven spatial distribution of different datasets leads to spatial gaps between the estimated 1D models. To overcome these gaps, we apply a Bayesian spatial fusion approach (Christakos, 1990; Serre et al., 1998) to construct continuous subsurface models. The estimated 1D models and their associated uncertainties are used to estimate the model parameters at the locations where no data have been measured. In the next section, the fusion workflow and its components are described. Then we illustrate the approach by a synthetic example. In Section 4 we apply the fusion algorithm to a field dataset with a specific focus on identification of soil contamination. 2. Bayesian fusion of multiple geophysical datasets Our fusion strategy is based on the capability of different geophysical methods to provide enough resolution to identify either subsurface material parameters or subsurface structure, or both. Combination of these two factors contributes to the overall information obtained from the geophysical investigations. In this paper we concentrate on three different geophysical methods: electrical resistivity tomography (ERT), electromagnetic induction (EMI), and GPR. Fig. 1 depicts the workflow of the fusion strategy where different geophysical datasets are prepared to enter a Markov chain Monte Carlo search-based inversion algorithm, and information content of the post-inversion results are quantified, thus allowing model appraisal and justification of the particular survey design. Spatial mapping of the target parameter is the final stage in the workflow, which leads to a coherent and smooth realisation of the target parameter. Although different geophysical methods are complementary, inconsistency in their acquisition geometries (i.e., 1D, 2D or 3D) makes the fusion process cumbersome. To overcome the incompatibility of the acquisition geometries and also to reduce computational cost of the 2D and 3D inversions (especially in a Bayesian framework), we re-cast the subsurface model to a set of locally layered 1D earthmodels. The assumption of a local 1D model is restricted to smooth geological structures and may not be appropriate for complex geologies such assumption may be employed to obtain initial subsurface models for 2D and 3D deterministic inversions. To estimate parameters of these 1D models, different geophysical datasets are treated in one of the following ways. 1. If the data are collected with 1D acquisition geometries, they may directly be imported into the inversion algorithm (e.g., EMI data). 2. If the data are collected with 2D or 3D acquisition geometries, they will be re-casted to equivalent 1D dataset (e.g., ERT data to VES data). 3. If the data are collected with 2D or 3D acquisition geometries and it is not possible to re-cast them to the equivalent 1D dataset, they will be interpreted separately and the result of the interpretation is used to build the prior model for the inversion. The co-located datasets are then jointly inverted to estimate the parameters of the 1D models. We jointly invert ERT and EMI datasets by using a joint likelihood function in the McMC search inversion. GPR data are processed independently to estimate the structural features, which are then used a priori to constrain the layers' thicknesses in the 1D models.

40

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Fig. 1. Geophysical data fusion workflow.

In contrast to the deterministic quasi-2D and 3D approaches mentioned earlier (e.g., Brodie and Sambridge, 2006; Viezzoli et al., 2008), we apply no smoothing or lateral constraints to the data. However, particularly in the 3D case where different data are collected at spatially dispersed locations, we apply a spatial fusion approach known as Bayesian maximum entropy (BME) (Christakos, 1990; Serre et al., 1998) to predict the model parameters at the locations where no data has been collected. The BME fusion approach uses the model parameters and their associated uncertainties estimated by McMC search-based inversion of the multiple datasets to produce a coherent 3D subsurface model. Herein, different components of the data fusion approach shown in Fig. 1 are explained. 2.1. Trans-dimensional model parameterisation Trans-dimensional inversion has become a popular inversion tool and has been successfully applied to a wide range of areas including inversion of DC resistivity data (Malinverno, 2002), multi-frequency EMI data (Minsley, 2011), receiver functions (Agostinetti and Malinverno,

2010) and travel-time tomography (Bodin et al., 2012). Here we follow the same approach that was proposed by Malinverno (2002) in the development of the inversion core of the fusion algorithm. The generic layered model that we use is shown in Fig. 2. The model is described by a model parameter vector m ¼ ðk; z; ρÞ;

ð1Þ

where k is the number of layers, z = (z1, z 2, z 3, …., z k-1) is a vector of depth to the layer interfaces and ρ = (log(ρ1), log(ρ2), log(ρ3), …., log(ρk)) is a vector of log-resistivities. These parameters, including the number of required layers, will be determined through the McMC sampling. Malinverno (2002) and Sambridge et al. (2006) show that because of the natural parsimony of the method, models with fewer layers will ultimately be dominant. We restrict the minimum zmin and maximum zmax depths based on the geometry of the measurement system and prior information. The number of layer interfaces is also restricted to fall between minimum kmin and

Fig. 2. Right: schematic representation of re-casted 2D model to a set of 1D earth-model. Left: 1D model parameterisation with k layers and corresponding assigned resistivity values. The minimum and maximum depths zmin and zmax, respectively, are specified by user based on acquisition geometry and prior information. The 1D model is adapted from Malinverno (2002) and Minsley (2011).

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

maximum kmax number of layers. The minimum thickness of layers is also restricted to: h min ¼

z max −z min : 2k max

ð2Þ

2.2. Bayesian formulation Non-linearity of the geophysical inverse problems and associated uncertainties in data and parameters can be analysed within a Bayesian framework. In this framework, the solution to an inverse problem is represented by a posterior (post-inversion) probability density function (pdf), p(m|d), over the model parameters m given data d (e.g., Tarantola, 2005). The posterior pdf is related to the so-called prior pdf, p(m|x), which represents pre-inversion knowledge about parameters m given prior information x, through: pðmjdÞ ¼

pðmjxÞpðdjm; xÞ pðdjxÞ;

ð3Þ

where p(d|m, x) is the likelihood function and the denominator in Eq. (3) can be regarded as a normalising constant because it is not a function of any model parameter (Sambridge et al., 2006). Prior knowledge about model parameters plays an important role in solving geophysical inverse problem. In the Bayesian formulation, the prior knowledge is presented through a multivariate pdf p(m|x). Malinverno (2002) and Minsley (2011) used multivariate Gaussian pdfs to describe the prior knowledge in similar inversion problems. However, we assume a uniform multivariate pdf within predefined ranges to describe prior information about log-resistivity, log ρ. Uniform prior pdf contains minimal information as compared to other types of prior models such as Gaussian distribution. Non-uniform prior pdfs priorities the model space to be inspected, however uniform prior pdf facilitates inspection of the model space without bias. The likelihood function describes how well the synthetic or modelled data, corresponding to each model m, match the recorded data, and here is defined to be a multivariate Gaussian in the difference between observed data diobs and synthetic data di(m): "  # 1 1 T −1 LðmÞ ¼ pðdjm; xÞ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  exp ∑ − 2 εj ∑j εj ;  j   ∏ 2π ∑j 

ð4Þ

j

where ∑ is the data error covariance matrix and ε is a vector representing the difference between observed data and synthetic data. Subscript j indicates different datasets (e.g., ERT and EMI). We assume here that data are independent and hence the covariance matrix is diagonal with elements σij for datum dij where subscript i = 1,2,… refers to the samples of jth dataset. For joint inversion of ERT and EMI datasets j = 2 and Eq. (4) is used as the joint likelihood. 2.3. Reversible jump Markov chain Monte Carlo search The goal of Bayesian approach is to characterise the posterior density p(m|d) in Eq. (3). It is reasonably well-known that finding an analytic solution to Eq. (3) is not easy, particularly in the trans-dimensional situation where the dimensions of the model space are large and variable. To overcome this, McMC sampling method has become a popular technique (e.g., Sivia, 1996). To obtain McMC samples of the trans-dimensional posterior pdf we use the so-called reversible jumb McMC (rj-McMC) (Green, 1995) which is an extension of the Metropolis-Hastings algorithm (Hastings, 1970; Metropolis et al., 1953). The sampling procedure contains two major stages of proposing a model in a probabilistic sense and then accepting or rejecting this proposal.

41

According to the Metropolis algorithm the probability of visiting a point in model space depends only on the current point and not on the previously visited points — this is the co-called Markov property. We sequentially update a set of estimates of the model parameters (e.g., log-resistivity and thickness described by vector m) using the following steps. • Initiate the McMC by calculating geophysical response d(mj) by forward modelling for an arbitrary model mj (i.e., a candidate number of layers and candidate depth to interfaces and log-resistivities). Then calculate the geophysical likelihood L(mj). • Define a new candidate parameter vector mj + 1 based on mj by randomly selecting candidate number of layers and then choosing one of the following jumps: (a) Birth with probability of 1/6: creates a new layer interface at a random depth. (b) Death with probability of 1/6: deletes a layer interface chosen randomly from the current set of interfaces. (c) Perturb with probability of 1/6: changes the depth of a layer interface chosen randomly from the current set of interfaces. (d) No change with probability of 3/6: leaves the layer interfaces unchanged. • The previous step will be completed by sampling the log-resistivity in each layer from the proposal pdf q(mj + 1|mj). • Calculate the corresponding geophysical response d(mj + 1) and likelihood L(mj + 1). • Use the Metropolis rule to accept or reject the new candidate model by calculating the acceptance probability:    3 L mjþ1 q mjþ1 =mj 5: α ¼ min41;   :  L mj q mj =mjþ1 2

If α is less than one, make a random decision to accept and set mj = mj + 1 or reject and set mj + 1 = mj. • Go to the first step • Repeat above steps until required number of samples in the set S = {m1, m2, …, mN} is obtained. The proposal pdfs for log-resistivity is multivariate Gaussian pdf as follows,  T     1 1 −1 exp − ; q ρjþ1 jk; z; ρj ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi − ρ ˜ C ρ − ρ ˜ ρ ρ jþ1   2 jþ1   ð2πÞkjþ1 C ρ  ð5Þ where Cρ is covariance of the proposed distribution. The mean vectors ρ ˜ is obtained from the log-resistivities in the current model. If the number of layers at step j + 1 did not change then ρ ˜ ¼ ρj . If the step j + 1 was a birth, the log-resistivities for the two layers that were generated by splitting a single layer are set the same as the original layer. In the case of death jump, the log-resistivity for the layer that is created by combination of two layers is the average of the corresponding values in the original layers. The general form of the probability of accepting a proposed model in trans-dimensional McMC may be written α ¼ min½1; prior ratio  liklihood ratio  proposal ratio  j Jj;

ð6Þ

where J is Jacobian that accounts for jumps between dimensions (Green, 1995). It has been shown that for the range of dimensional changes that we are dealing with the Jacobian term is equal to one (Agostinetti and Malinverno, 2010). Since we are using uniform

42

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

prior distributions, the prior ratio is also equal to one. Hence, our expression of α becomes, α ¼ min½1; liklihood ratio  proposal ratio:

ð7Þ

The McMC algorithm will start sampling the posterior after a ‘burn-in’ period. The problems associated with choosing the burn-in period and the maximum number of samples in McMC algorithms are on-going subjects of research (Cowles and Carlin, 1996; Gelman and Rubin, 1992; Gilks et al., 1996). Practically, the burn-in period may last until the misfit between data and model response falls below the data error range and the McMC may continue until the characteristics of the posterior pdf show no significant change (Malinverno, 2002; Minsley, 2011). The set S contains a set of samples that approximate samples of the posterior pdf. By calculating the sample density in S, we obtain an estimate of the posterior pdf p(m|d) for any new measured data dobs. The speed of sampling the posterior in McMC strongly depends on the choice of proposal distribution. In order to improve the efficiency of sampling the posterior, Gelman et al. (1996) recommended keeping the acceptance rate τ between 0.23 and 0.44. If τ > 0.44, the proposal variance is increased (e.g., by 10%) in order to decrease the acceptance rate and if τ b 0.23, the proposal variance is decreased (e.g., by 10%). Gelman et al. (1995) use this adaptive approach to tune the McMC algorithm and then run the algorithm with a constant proposal distribution to sample the posterior distribution. Roberts et al. (1997) also show that the acceptance rate optimizing the efficiency of the McMC approach is 0.234. Hereon we use the Gelman et al. (1996) approach by calculating τ for every 50 iterations throughout the Markov chain and adjusting the proposal variance. This adaptive approach may not be the most efficient method to improve convergence of the McMC towards correct posterior distribution. In this paper we do not intend to optimize the computational efficiency of the McMC sampler, although we appreciate that methods have been developed to improve the acceptance rates within the Metropolis algorithm beyond 0.23 (Bedard, 2008) or for adaptively updating the scale and orientation of the proposal distribution (Vrugt et al., 2009). There is usually a burn-in period at the beginning of the Markov chain before it starts sampling the posterior distribution. We specify the end of the burn-in period when a model in the Markov chain fits the data within an expected error tolerance (i.e., the predefined data error). The length of this period depends on the assumed data error and it may last for several hundreds or even few thousands of model proposals if the assumed data errors are very small. Assigning a proper length to the McMC sampling is particularly difficult for trans-dimensional algorithms. We run the Markov chain for tens of thousands of steps until the variance of the posterior distribution of the top-layer's log-resistivity does not change significantly. Although multiple Markov chains may be run in parallel to help ensure that the whole model space is inspected, we assume that in our case a single chain is converging towards the global solution. This assumption is supported by running a number of McMC for a synthetic VES example with different starting models. 2.4. Information quantification The state of knowledge about any random or uncertain parameter that is described by a pdf may be quantified by using an information measure (Lindley, 1956). Shannon (1948) introduced a unique, linear measure of information represented by any pdf f(x) by: I½f ðxÞ ¼ −Ent ðxÞ þ c ¼ ∫ f ðxÞ log ½ f ðxÞ dx þ c:

ð8Þ

x

Here I is the information measure as defined by Shannon (1948), Ent is the entropy function defined in the right hand side, f(x) is the

pdf of the random variable x, and c is a constant. Shannon’s entropy has been widely used in experimental design theory and applications (e.g., Curtis, 2004a,b; Guest and Curtis, 2009, 2010; JafarGandomi and Curtis, 2012; Krause et al., 2008; Sebastiani and Wynn, 2000; van den Berg et al., 2003; Zidek et al., 2000). The unit of information depends on the base of the logarithm in Eq. (8). The corresponding units for logarithm bases 10, Euler’s number e and 2 are dit, nat and bit, respectively. From here on, we use base-e logarithm for estimation of information. Generally speaking, information about a parameter is maximized when its pdf constrains its possible ranges of values most tightly (Lindley, 1956). It is useful to estimate Shannon's information of a pdf such as f(x) with respect to another pdf such as g(x) by using a relative information measure (Cover and Thomas, 2006),  f ðxÞ dx þ c: Irel ½ f ðxÞjgðxÞ ¼ −Ent rel ðxÞ þ c ¼ ∫ f ðxÞ log g ðxÞ x

ð9Þ

The calculation of Shannon's entropy is a significant issue from a computational point of view. Ahmed and Gokhale (1989) showed that the Shannon's entropy of a multivariate Gaussian distribution is Ent = log[(2π e)N|Σ|]/2 where N indicates number of variables. Shannon's entropy of a uniform distribution is log |b − a| where a and b are lower and upper bounds of the distribution. Cover and Thomas (2006) give an analytical solution for Eq. (8) for a particular set of pdfs. However, in many real world problems it is not possible to represent uncertainties with any single analytical pdf. An approximation that can be demonstrated to converge to the true entropy relies on random sampling methods such as Monte Carlo used in this study. To calculate the information we use the normalised histogram (with unit volume) of the estimated models, which is considered as a discretised numerical approximation to the posterior pdf. Entropy has several advantages over commonly used statistical measures to quantify pdfs. For example, while the variance depends just on the second moment of pdf, entropy has the advantage of depending on many more parameters which are embedded in the pdf's shape. Krykacz-Haussmann (2001) shows that the conditional variance (variance of a conditional pdf) may be larger than one, while this is not the case for entropy. 2.5. Spatial fusion (mapping) For a particular site, when spatial heterogeneity of the target soil parameter (in our case, electrical resistivity) is relatively low and sampling/measurement locations are regularly distributed over the site, basic interpolation algorithms such as Kriging (e.g., Cressie, 1985) might be pragmatically acceptable for spatial prediction of the target parameter at un-sampled locations. However, quite often at least one of the above assumptions are disturbed, hence a different approach must be applied. Statistical spatial data fusion is becoming a popular tool for spatial prediction based on dispersed observations with different levels of uncertainties (e.g., Bogaert and Fasbender, 2007; Clemen and Winkler, 1993; Douaik et al., 2004; Fasbendar et al., 2009; Jouini and Clemen, 1996). We adapt a Bayesian spatial prediction approach similar to that employed by Douaik et al. (2004) which is based on the Bayesian maximum entropy (BME) (Christakos, 1990; Serre et al., 1998). BME is a spatial (and temporal) prediction method that can incorporate a wide variety of datasets with different uncertainties. Using BME, the target subsurface parameter is considered as a random variable. A specific realization of the random variable is χmap = (χhard,χsoft,χk), which consist of three components; the hard data (accurate) χhard, the soft data (vague) χsoft, and the data to be estimated χk at locations Phard, Psoft, and Pk, respectively. Note that χhard and χsoft are known and χk is unknown and will be estimated through BME. The estimation locations Pk may form a regular grid to be used to produce a map of χk. In order to maximize the

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

information about the target random variable over a certain area, three stages are introduced: 1. At the first or prior stage, an initial probability distribution will be generated across space on the general knowledge base (G), which includes physical laws as well as other sources of information, as available. Such information may be included in G as mean trend and the covariance function for the target parameter within a specific region. 2. At the second stage (meta-prior), the site-specific knowledge available is organized into hard (accurate) and soft (vague) data and is denoted by S. 3. At the third or integration (posterior) stage, the site-specific knowledge of stage 2 (S) is assimilated with G. The final solution contains the complete probability law at each estimation point. The mathematical formulation of these three stages is as follows. Consider a geophysical parameter such as electrical resistivity X(s) obtained form joint McMC inversion of geophysical data for a particular depth over the target site. This parameter is constrained by the general knowledge G as     G smap ¼ g^ α xmap ;

α ¼ 1; 2; …; NC

      g^ α xmap ¼ ∫g α χ map f G χ map dχ map ;

ð10Þ ð11Þ

where the vector xmap indicates the random variables associated with X(s) at a set smap of mapping points, gα(χmap) are a set of functions of χmap whose expectation is known from G (Serre and Christakos, 1999) and facilitates incorporation of the general knowledge. In our case, the general knowledge G consist of the mean mX ðsÞ ¼ E½X ðsÞ;

ð12Þ

and covariance function C ðs; s′ Þ ¼ Ef½X ðsÞ−mX ðsÞ½X ðs′ Þ−mX ðs′ Þg;

ð13Þ

where E[.] indicates statistical expectation. fG(χmap) is the G-based multivariate pdf representing the prior knowledge, which has to be found. The general knowledge G is used to maximize the expected amount of information, which is represented by Shannon’s entropy Ent:   h  i Ent ¼ −∫f G χ map log f G χ map dχ map ;

ð14Þ

Eq. (14) may be maximised by introducing the Lagrange multiplier μ α: h  i   h  i L f G χ map ¼ −∫f G χ map log f G χ map dχ map h     h  ii −∑μ α ∫g α χ map f G χ map dχ map −E g α χ map ; ð15Þ the values of μα are obtained by solving Eqs. (14) and (15). Solving the system of Eq. (14) leads to   h  i f G χ map ¼ Z exp ∑μ α g α χ map ;

ð16Þ

where Z is a normalisation constant (see Christakos, 2000, p. 74–76). Knowing μα and hence the G-based prior pdf (Eq. (16)) means that we have completed the first stage. At the second stage the site specific knowledge S is organised as hard and soft data. The hard data present very high reliability and accuracy and the soft data indicate lower accuracy or precision, and are usually associated with a certain level of uncertainty.

43

At the third stage, the general knowledge G (i.e. prior pdf) is updated with the site-specific knowledge S to reach to the total knowledge K = G ∪ S, which is represented as the BME pdf given by:   h  i K −1 f K χ ¼ ðAZ Þ ∫ exp ∑μ α g α χ map dχ soft ;

ð17Þ

where A is a normalisation coefficient. The BME pdf is in general non-Gaussian and incorporates spatial dependency and accuracy of the data. The crucial advantage of the BME prediction over Kriging is distinguishing between hard and soft data. The soft data may be in the form of pdfs or intervals. Serre (1999) provide the following solution for the pdf-type (Eq. (14)) and the interval-type (Eq. (15)) soft data, respectively,        ∫f G χ k ; χ hard ; χ soft f S χ soft dχ soft     ¼ ; f K χ k jχ hard ; f S χ soft ∫f G χ hard ; χ soft f S χ soft dχ soft

ð18Þ

     ∫βα f G χ k ; χ hard ; χ soft dχ soft   f K χ k jχ hard ; f S χ soft ; ¼ β ∫α f G χ hard ; χ soft dχ soft

ð19Þ

where α and β are lower and upper bounds of the intervals. Note that for all BME computations we use the BMElib program package which is a MATLAB numerical toolbox for spatio-temporal geostatistical analysis based on the BME approach. Fig. 3 shows an example of application of BME to recover the true synthetic model (Fig. 3a) from 60 sampled data. The data consists of 10 hard data samples and 50 soft data samples. The hard data (asterisks) are represented by the exact value of the true model at the sampling locations. Soft data (circles) are also taken from the true model and are represented by ±10% of the true values at sampling locations (interval-type). Fig. 3b and c depict the predicted models using BME method and Kriging method, respectively. The BME estimate is more successful in recovering the true model than the Kriging estimate. In particular, boundaries of the estimated model using BME (Fig. 3b) are much sharper than those for the estimated model using Kriging (Fig. 3c). 3. Synthetic case study DC resistivity, EMI and GPR methods are recognized to be valuable methods for characterization of hydrologic systems due to their high sensitivity to the water content, fluid salinity and textural properties of shallow soils. EMI instruments have been successfully used for soil characterization and estimation of hydraulic parameters (e.g., Brosten et al., 2011; Cameron et al., 1981; Hezarjaribi and Sourell, 2007; Triantafilis and Monteiro Santos, 2011). Resistivity methods, such as ERT and VES, have long been used in water resource exploration and hydrologic system characterization (e.g., Crook et al., 2008; Kosinski and Kelly, 1981). In recent years, GPR has also become a popular tool for near-surface site characterization (e.g., Dubreuil-Boisclair et al., 2011; Huisman et al., 2003; Meles et al., 2012). Here, we focus on these popular methods to demonstrate performance of the fusion approach with a synthetic example. We start with a simple example for a 1D earth model. The synthetic model contains two layers with 2 m and 10 m thicknesses laying on a half-space. The resistivity values of the layers and the half-space are 100, 1000 and 100 ohm m, respectively. To produce synthetic data, we forward model VES and EMI data on the surface using the approaches developed by Ingeman-Nielsen and Baumgartner (2006) and McNeill (1980), respectively. A review of these forward modelling approaches is provided in Appendices A and B. We add 2% random noise with Gaussian distribution to the synthetic VES and EMI data. We assume that GPR data are also collected on the surface and are able to provide us with the depth to the top interface with 20% uncertainty. To generate the VES data, we assume a dipole–dipole electrode

44

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Fig. 3. Comparison of capability of BME method and Kriging method to recover (a) the synthetic model. Open circles and asterisks indicate soft and hard data, respectively. (b) Estimated model by using BME method and (c) Kriging method.

configuration with dipole size of 10 m. The synthetic sounding consists of 20 apparent resistivity values at 20 dipole spacings from 2 m to 97 m. We generate a three-sample EMI dataset corresponding to EMI measurements in a vertical dipole mode with coil separations of 1 m, 3.66 m and 7 m. These datasets are inverted individually and jointly using the McMC search-based inversion algorithm with 100,000 samples, which is required to obtain a stable posterior distribution, and burn-in period of 2000 samples. We select one sample out of every 50 samples to construct the posterior distribution. The prior model for the log-resistivity of layers is uniform between 1 and 4.3 corresponding to 10 ohm m and 20,000 ohm m, respectively. The thickness of the top layer is sampled from uniform distribution centred on the estimated depth from GPR data bounded with ±10% of the centre value (2 m). The minimum and maximum allowed depth for interfaces are 0.25 m and 25 m, respectively. Fig. 4 depicts the impact of fusion of different datasets on the estimated subsurface model. Fig. 4a shows the posterior distribution of the estimated resistivity models at different depths obtained from inversion of VES data. To construct this figure the selected models from the McMC inversion are discretised with 10 cm depth-intervals and then at each depth the normalised histogram of the models is calculated to produce density distributions. The colour scale indicates density of the estimated resistivity values. Solid and dashed red lines indicate median and 95% confidence interval of the estimated models and the solid white line indicates the true model. Fig. 4b shows the posterior distribution of the estimated models from joint inversion of VES and EMI datasets. This figure shows that joint inversion of the VES and EMI data greatly improves the estimated subsurface model. Note that the improvement is not limited only to the shallow part of the model. Additional information from the EMI dataset improves also the deeper parts of the model. Fig. 4c indicates impact

of fixing the top layer thickness (and the minimum depth zmin) by bringing-in the information from GPR, which significantly improves the estimated model. We calculate relative information of the estimated posterior distributions with respect to the prior distribution (Eq. (9)) from different inversions (Fig. 4d). This figure implies improvement of the relative information content of the estimated posterior distributions by including additional datasets into the fusion algorithm. 4. Field case study We now apply the fusion tool to a collection of datasets obtained from a field site located close to Trecate (Novara – NW Italy). At this site multiple surface and borehole-based geophysical measurements have been conducted to characterise the near surface hydrogeology and to identify the contamination plume caused by a nearby oilwell blowout in 1994. 4.1. Site description The Trecate site is mainly used for agricultural purpose and the near surface layers are characterized by a sequence of poorly sorted silty sands and gravels in extensive lenses, typical of braided river sediments. An artificial layer of clayey-silty material, 1–2 m thick originally placed as a liner for rice paddies overlies most of the site. This succession can be clearly seen on the lithology logs and electrical conductivity logs (EC-logs) obtained by direct-push logging tool. At this site the total petroleum hydrocarbon (TPH) in the soil has been measured at different depths with direct sampling (Eni, 2008). Fig. 5 depicts interpolated map of the depth-averaged TPH distribution over the site.

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

45

Fig. 4. Inversion results of the synthetic data. (a) Posterior distribution of individual McMC inversion of VES. The colour scale indicates density of the estimated models, which are calculated by discretizing the estimated models with 10 cm depth-intervals and calculating the normalized histogram at each depth. (b) Posterior distributions of the joint inversion of VES and EMI data. (c) Posterior distributions of the joint inversion of VES and EMI data constrained with information from interpretation of the GPR data. Dashed lines indicate 98% confidence interval of the posterior distributions and solid line shows the posterior median. (d) Depth-dependent variation of relative information estimated from posterior distributions of individual and joint inversion of the synthetic VES, EMI and GPR data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Interpolated map of the average TPH values obtained from direct sampling. Rectangles indicate location of the boreholes where EC-logs have been collected and cross symbols show the sampling locations for TPH measurement.

46

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Fig. 6a–g show the density distribution of the logarithm of resistivity at different depth intervals extracted from the EC-logs at the contaminated and non-contaminated regions of the site. The locations of the boreholes are depicted in Fig. 5. These plots indicate that, particularly at the depth interval of 2 m–10 m, there is a significant shift towards lower resistivity values at contaminated regions. At the depth interval 10 m–12 m, which coincides with the water table depth, both probability distributions are shifted towards lower resistivity values. Below the water table there is no clear distinction between the resistivity values at the contaminated and non-contaminated regions. In this case, a decrease of soil resistivity appears mainly due to the increased microbial activities in the vadose zone in presence of the mature plume of petroleum hydrocarbon, which is well-evident in the literature (e.g., Atekwana and Atekwana, 2010). The potential high sensitivity of the soil resistivity to TPH makes it an appropriate parameter for further site characterization using geophysical methods. 4.2. Geophysical data At the Trecate site, 10 ERT and 16 GPR transects, 21676 EMI measurements using an EMI instrument with coil-separation of 3.66 m (EMIDeep) and 1061 EM measurements using an EMI instrument with coil-separation of 1 m (EMIShallow) were collected. The contour-maps

in Fig. 7a and b show the interpolated EMIDeep and EMIShallow datasets, respectively. Locations of the ERT and GPR transects are shown in Fig. 7a and b, respectively, with dashed lines. ERT data are retrieved by using dipole-dipole electrode configurations with electrode spacings of 1.5 m, 2.5 m, and 3 m for different transects. The dipole-size is equal to 4 electrode spacings. GPR data are collected with a common-offset acquisition and constant transmitter-receiver offset of 1 m. Fig. 7c shows an example of the recorded GPR data along Transect B. A direct-push tool was also used to drill 16 boreholes for collecting EC-logs. Locations of these boreholes are shown with crossedrectangles in Fig. 7a. A range of geophysical methods including ERT, EMI, and GPR were selected to characterize the site based on their sensitivity to the soil resistivity and/or subsurface structural features. The spatial coverage and resolution of the applied methods are very different, which makes the fusion process cumbersome. Such situations usually occur due to acquisition inconsistency of the methods, practical limitations and physical barriers in the field. Prior knowledge and targeting specific regions of the field may also lead to non-uniform spatial coverage. Hereon, strategies for fusion of different datasets in two- and three-dimensional spaces are developed. The EMI data were collected in a point-wise manner and can be readily imported into the fusion algorithm. Following the fusion

Fig. 6. (a)–(g) Probability distribution of the soil resistivity (shown in log10 (resistivity, in ohm m)) obtained from EC-logs at contaminated (solid line) and non-contaminated (dashed line) regions of the site.

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

47

Fig. 7. (a) EMIDeep data, (b) EMIShallow data and (c) an example of GPR transects collected at the Trecate site. Location of ERT and GPR transects are shown with dashed lines and red circles in (a) and (b), respectively. Crossed rectangles in (a) indicate locations of the boreholes drilled by direct-push tool where EC-logs were obtained.

algorithm in Section 2, we re-cast the 2D ERT data to a set of 1D VESs by re-arranging the measured apparent resistivities to be inverted for the locally 1D earth-model at each sounding location. The 2D GPR data are processed to extract the depth to the layer interfaces at each sounding location. 4.3. Data fusion along cross-sections (quasi-2D) At the Trecate site, while the EMI data were collected with high spatial coverage, the ERT and GPR data were recorded along 2D transects with un-even distribution over the site. Although, the measurement frequencies of ERT and EMI methods are different, our modelling based on that of Pride (1994) shows that the frequency dependent variation of electrical resistivity within our measured range of frequencies is negligible. The two EMI dataset are collected with different instrument heights, which are considered in the inversion. Fig. 8 depicts fusion of EMIDeep, GPR and ERT (VESs) datasets along Transect B (the location of Transect B is shown in Fig. 7a). Approximate spatial locations of the measured apparent resistivities along Transect B used to construct VESs are shown in Fig. 8a. Fig. 8b shows interpolated EMIDeep data and estimated depth to the first interface extracted from the GPR data at each VES location. Importing these datasets into the fusion tool provides the subsurface log-resistivity image (Fig. 8c), associated uncertainties (standard deviation) (Fig. 8d) and relative information (with respect to the prior distribution) images (Fig. 8e).We use a

uniform prior distribution for the log-resistivity between 1 and 4.3 corresponding to 10 ohm m and 20,000 ohm m, respectively. We assume Gaussian data error with standard deviation of 40 ohm m for the VES data (apparent resistivity values) and 0.3 mS/m for the EMI data (apparent conductivity). Note that using the same data error for both low and high apparent resistivity values may affect the estimated posterior distributions. We let the number of layers vary between 2 and 30. The maximum depths for the interfaces is variable along the transect and is dictated by the maximum electrode-spacing for each sounding location. The thickness of the top layer is constrained with the GPR data and is also considered as the minimum allowable depth for the interfaces. Note that in this case we did not apply the spatial mapping (BME) part of the fusion algorithm (Fig. 1) along the cross-section because of the close proximity of the extracted VESs. EC-logs are available at three locations on Transect B and indicate a relatively good agreement with the inversion results, however, at the deeper interface at the middle of the cross-section agreement with the EC-logs is weak. The EC-logs and the lithological logs obtained from the boreholes indicate a general three-layer geological model along the cross-section and also for the rest of the site. Therefore, two main resistivity contrasts are expected to be seen on the estimated resistivity image. Reliability of the estimated resistivity values along the cross-section at different depths is clearly reflected in the estimated information (Fig. 8e). Generally, along the cross-section relative information decreases with depth. However,

48

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Fig. 8. Fusion of EMIDeep, GPR and ERT dataset along Transect B. (a) Approximate spatial location of measured resistance along the ERT transect, (b) Interpolated EM measurements and estimated thickness of the top-layer extracted from interpreted GPR transect, (c) mean of the estimated log-resistivity (ohm.m) value using quasi-2D McMC inversion, and corresponding (d) standard deviation and (e) relative information for estimated resistivity values. EC-logs obtained from direct-push tool are also shown for comparison.

because at the edges of the cross-section the maximum allowable depth is tightly constrained, this constraint increases the level of relative information. Further inspection of the estimated relative information indicates higher values of information for the middle layer. This is because of the proximity of the resistivity values of the middle layer to the upper bound of the prior distribution. JafarGandomi and Curtis (2012) show that bounds of the prior distribution act as hard-boundary conditions and increase the level of obtainable information. Variation of standard deviation of the estimated models along the cross-section indicates a general trend of increasing uncertainty with depth and also higher uncertainty around the major resistivity contrasts. The shallower contrast is obviously well-resolved due to the constraint posed by the GPR data. The deeper contrast is associated with higher uncertainty. This uncertainty may be reduced by incorporation of the EC-logs and the lithological logs in the fusion, which were not include in the current inversion. As mentioned earlier, the assumption of locally layered earth for inversion of VES data will be violated in the case of complex geological structures and particularly large offsets. However this does not seem to significantly affect the inversion results at the Trecate site due to its relatively simple near-surface structure. Fig. 9 shows comparison of the McMC inversion and deterministic inversion results for a selected VES data from Transect B (VES-B30 in Fig. 8a). This figure indicates a significant agreement between the deterministic model and the most probable model from McMC inversion.

4.4. Data fusion in three-dimensions In order to integrate the dispersed datasets in 3D, we use a quasi-3D approach similar to the quasi-2D approach in Section 4.3. We introduce a spatial grid with 10x10m grid-cells and average the collected data for each method within the grid-cells. We assume that the subsurface model within each grid-cell is vertically heterogeneous with no lateral variations. All the ERT transects are re-casted to VESs and the GPR transects are interpreted to extract the top layer thickness. The VES data within each grid-cell are averaged for each individual dipole-spacing to construct a single averaged-VES for each cell. Each grid-cell may contain data from one, two, three or none of the datasets. A large number of the grid-cells contain only EMIDeep data, which means that no information about vertical variation of the soil resistivity can be obtained at these grid-cells. However, we use the prior knowledge about the approximate thickness of the top soil layer to constrain the models. We use the estimated thickness of the top layer at each grid-cells along the GPR transects to map the top-layer thickness over the entire site using the BME approach. We consider the GPR dataset as soft data with 10% uncertainty in the thicknesses extracted from interpretation of the GPR transects. The covariance model of these data is estimated by fitting an exponential function to the extracted thickness. BME provides us with the estimated thickness and corresponding uncertainty at each grid-cell. Fig. 10a and b show variation of the thickness of the top-soil layer and associated

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

49

Fig. 9. An example of the inversion results for a VES on Transect B. The location of this VES is labelled with VES-B30 in Fig. 8. (a) Apparent resistivity data and the best fit curve (solid red line) obtained from deterministic inversion of the data. Colour scale indicates density (from normalised histogram) of the apparent resistivity curves predicted by the layered media sampled by the McMC algorithm. (b) Posterior distribution of individual McMC inversion of VES-B30 data. The colour scale indicates density of the estimated models. Dashed lines indicate 98% confidence interval of the posterior distributions and the solid red line shows the deterministic inversion result. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

uncertainty (variance) obtained from BME fusion of the estimated thickness along the GPR transects. We ran the McMC search-based inversion at each grid-cell. The prior distribution for the resistivity was uniform between 1 and 4.3 corresponding to 10 ohm m and 20,000 ohm m, respectively, and the thickness of the top layer is constrained by the GPR information (Fig. 10). The maximum depth is constrained by the maximum dipole-spacing for the soundings at each grid-cell. For those gridcells that only EMI data exists, the maximum depth is 4.5 m. the number of interfaces may vary between 2 and 30. The number of samples for the McMC inversion is 100,000 with a burn-in period of 2000 samples; we selected one sample out of every 50 samples to estimate the posterior pdf. Fig. 11 shows maps of the averaged estimated relative information for the resistivity in three depth-intervals of 0–1 m, 3.5–8.5 m and > 12 m at those grid-cells where McMC inversion is conducted. In the south-east of the site where both EMIDeep and EMIShallow datasets are available the information is particularly greater for the top

depth-interval (0–1 m), as expected intuitively. Including the ERT dataset significantly improves the information, particularly where no EMIShallow data are available. Because of the limited investigation depth of the EMI measurements, the information at the deeper intervals is available only for the grid-cells with ERT measurements. For most of the grid-cells where only EMIDeep data are available the estimated information for the top-interval is correlated with the thickness of the top layer. This is explained by the cumulative depth-sensitivity (Eq. (A2)) of the EMI measurement (e.g., McNeill, 1980), for which 70% of the measured EMI response corresponds to the top 6 m of the soil. Fig. 11d shows three examples of the estimated posterior distributions for the three depth-intervals. These distributions are obtained by averaging the posterior distributions over each depth-interval at three different locations labelled with 1, 2, and 3 in Fig. 11a, b, and c, respectively. This figure indicates the non-Gaussian characteristics of the posterior distributions, which are taken into account for spatial mapping of the estimated resistivity values.

Fig. 10. (a) Map of estimated thickness of the top layer obtained from BME fusion of extracted thicknesses along the GPR transects and (b) corresponding estimated variance.

50

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

The results in Fig. 11 imply that fusion of the EMI and GPR datasets provides valuable information for the top-interval, however, for the middle-interval the ERT dataset is the main contributor to the information. Fig. 11 provides valuable information from a survey design point of view, which may be used to develop further geophysical surveys at the site. For example, while ERT method provides a relatively higher level of information, its low spatial coverage may be compensated by using the EMI method. Although, the EMI measurements with single depthsensitivity provide great coverage, their sensitivity to depth variation of resistivity may be improved by including additional EMI measurements with different depth sensitivity and also by integrating them with the GPR dataset. Identifying the optimal combination of the above mentioned methods depends on their availability and cost combined with their capability to provide subsurface information. After estimating the 1D subsurface models at all grid-cells with data, we create depth slices over the site to build the 3D model. Fig. 12 shows three depth-slices (maps) from the estimated 3D model. To produce these maps, estimated resistivity models in the form of pdfs (i.e., normalized histograms) at each depth-interval are fed into the BME fusion algorithm to predict the model parameters on a regular grid with 10 m grid spacing. Using the estimated resistivity values in the pdf forms guarantees incorporation of non-Gaussian characteristics of the

posterior distributions into the spatial mapping. We obtain the covariance model for each depth interval (which is needed for BME estimation) from the EC-log data at the corresponding depth intervals. We assume that such covariance model represents the lateral covariance of the geological structure. The best-fit exponential covariance models of the top, middle and bottom intervals are characterised with sill values of 0.1, 0.3 and 0.3 (in log-resistivity units) and the spatial ranges of 150 m, 50 m, 30 m, respectively. All estimated models are considered as pdf-type soft data. We use up to five data within a 30 m radius from each grid-node to estimate the corresponding log-resistivity and associated uncertainty. The left and right panels in Fig. 12 show maps of the estimated log-resistivities and uncertainties (variance), respectively, at the three depth-intervals of 0–1 m, 3.5–8.5 m and >12 m over the entire site. The uncertainty in the estimated resistivity (Fig. 12 right panels) is directly related to the distance from the measurement locations and the uncertainty in the measurements, which is reflected in the estimated posterior pdfs and covariance model. The log-resistivity maps are best interpreted along with the uncertainty (variance) maps, particularly when decisions have to be made for further investigations and to design subsequent surveys. To identify correlation between the estimated resistivity values and the contamination, we compare the estimated resistivity values

Fig. 11. Depth-average of the estimated relative information from fusion of EM, GPR and ERT datasets for the resistivity of (a) top depth-interval (0–1 m), (b) middle depth-interval (4.5–10 m) and (c) deep depth-interval (>10 m).

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

51

Fig. 12. Maps of estimated log-resistivity for (a) the top depth-interval (0–1 m), (c) middle depth-interval and (e) deep depth-interval and their corresponding estimated variances, (b), (d), and (f), respectively. These maps are derived from BME spatial fusion of the estimated 1D models. All resistivity scales are in log10 (resistivity, in ohm m).

at the known contaminated grid-cells with those at the known non-contaminated grid-cells for the top two depth-intervals (0–1 m and 3.5–8.5 m). The contaminated and non-contaminated grid-cells are identified by using the TPH distribution map (Fig. 5). Fig. 13 shows the corresponding density distributions (normalized histograms). For the second depth-intervals the density distribution for the contaminated grid-cells shows a clear shift towards lower resistivity values. This result is in agreement with that obtained independently from the EC-logs shown in Fig. 6, which were not used for inversion. Comparison of the results shown in Figs. 6 and 13 indicates

successful application of the fusion tool at the site for mapping soil resistivity using multiple geophysical datasets. 5. Computational constraints Although the proposed fusion algorithm is a significant step towards feasibility of 3D characterisation of subsurface with large multiple datasets, the locally-layered assumption for building the 3D model may not be appropriate in highly complex geological models. However, even in such cases the Bayesian fusion algorithm is greatly helpful to

52

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Fig. 13. Density distribution of the estimated log-resistivity (ohm m) values at the contaminated and non-contaminated grid-cells for (a) the depth-interval (0–1 m), (b) the middle depth-interval (3.5–8.5 m). The contaminated and non-contaminated grid-cells are separated by using the boundaries of contaminated plumes shown in Fig. 5.

obtain an initial 3D model which may be used for later focused 3D inversions. Furthermore, the spatial distribution of estimated parameter uncertainty is a valuable outcome of the fusion algorithm that may be used to design further experiments and to identify the appropriate methods, which help the most to reduce subsurface uncertainties. An issue of the McMC methods is their computing time requirement. Usually a very large number of samples are required to obtain a posterior distribution that best represents the quality of the estimated model parameters. For example, in our case running the McMC search-based inversion for a single VES with 15 data samples and 50,000 iterations takes 10 min on a computer with a dual-core 3 gigahertz processor and 4 gigabyte memory. Choosing a proper proposal distribution can significantly increase efficiency of the posterior construction. Gelman et al. (1995) proposed adjusting the variance of the proposal distribution frequently along the Marko chain — the method that we use in this paper, to improve the convergence and efficiency of the search. Tierney and Mira (1999) proposed the delayed rejection scheme to improve convergence rate of the Markov chains. This approach was used by Green and Mira (2001) for trans-dimensional moves within McMC and by Bodin and Sambridge (2009) in seismic tomography. According to the delayed rejection scheme, once a rejection decision is made, instead of moving to the next iteration and remaining in the same position in the parameter space, the variance of the proposal is decreased and a second move is proposed. Improved computational performance (hardware) coupled with more efficient computational algorithms, such as message passing interface (MPI), will clearly reduce computational time. Since the fusion algorithm consists of a large number of independent individual 1D inversions, which may be run in parallel, the overall computation cost will significantly be reduced when a large number of computer cores are used simultaneously. In spite of the computational issues, McMC is becoming a popular method for Bayesian inversion of geophysical data. Note, however, the posterior probabilities, upon which the data fusion methodology relies, are strongly dependent on the validity of the local 1D approximation and the specific data error structure assumed in the McMC procedure. An alternative approach to implement the fusion algorithm may be developed by replacing the “stochastic joint inversion-McMC” part of the proposed fusion algorithm (Fig. 1) with a set of deterministic 2D inversion of ERT, EMI and GPR dataset. The estimated resistivity values then may be input into the BME fusion procedure. Although relative accuracy of these two approaches may be investigated, nevertheless the fact that the McMC inversion provides non-Gaussian posterior distributions, which may be used for BME fusion, is a key advantage of the proposed approach.

6. Conclusions We propose an approach for fusion of multiple geophysical datasets within a Bayesian framework. Our approach is developed within the framework of characterising near surface soils. In the fusion workflow, different geophysical datasets are fed into a trans-dimensional Markov chain Monte Carlo — based joint inversion algorithm and the information content of the post-inversion results are quantified with Shannon's information measure. To improve efficiency and practicality of the algorithm, we re-cast the model space to a set of closely-spaced locally 1D layered models, which affects the estimated posterior probabilities at sites with complex geology. We use a Bayesian maximum entropy (BME) approach for spatial mapping of the target parameter. The fusion approach is capable of integrating geophysical methods with different spatial coverage and resolution by using quasi-2D and quasi-3D approaches to represent subsurface models. Computation of information contents permits an assessment of data worth for each geophysical modality: an aspect that has been hitherto overlooked in the majority of near surface geophysical investigations. Application of the fusion tool to a hydrocarbon contaminated site indicate that the proposed fusion strategy has potential for not only in enhancing the subsurface information but also as a survey design tool to identify the appropriate combination of the geophysical tools and thus show whether application of an individual method for further investigation of a specific site is beneficial. In this case, although the electrical resistivity soundings provide a high level of information about the subsurface, its low spatial coverage is compensated by using electromagnetic induction method. Based on the results of the data fusion tool, and independently from EC-logs, electrical resistivity appears to be a useful indicator of soil contamination in the vadose zone at the study site. The fused geophysical model may then prove useful for focussing, and monitoring, targeted remediation, in areas where direct observations are not available. We anticipate further application of this practical fusion tool, particularly in respect to the optimisation of geophysical surveys for site characterisation. Acknowledgments This research was made possible by funding from the EU FP7 collaborative project ModelPROBE ‘Model Driven Soil Probing, Site Assessment and Evaluation’. The data used in the project are a result of significant input from colleagues at the University of Padova (Italy) and the University of Bonn (Germany). In particular, the authors acknowledge contributions from Adrian Flores-Orozco, Giorgio Cassiani

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

and Andreas Kemna. We also wish to thank editor Prof. Klaus Holliger and two anonymous reviewers for valuable comments. Appendix A The electromagnetic induction (EMI) instruments induce a primary magnetic field from a transmitter coil into the ground, which leads to a secondary magnetic field recorded at the receiver coil. Since at low induction numbers (LIN) condition, the magnetic couplings between induced eddy currents are negligible, the received secondary magnetic field is the sum of all ground induced ground current loops (McNeill, 1980). Kaufman and Keller (1983) indicate that in this condition transmitter-receiver separation is the sole factor affecting the effective depth of exploration. Assuming a horizontally layered earth-model, McNeill (1980) developed a cumulative sensitivity function (CS) representing relative contribution of the subsurface materials up to a depth d below the sensor to the secondary magnetic field: CSHMD ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4z2 þ 1−2z;

ðA1Þ

for a horizontal magnetic dipole (HMD) and, 1 CSVMD ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 4z2 þ 1 for a vertical magnetic dipole (VMD), where z indicates normalized depth (actual depth divided by transmitter-receiver distance). Based on the cumulative sensitivity relations, the expected response of an N-layer earth-model in terms of apparent conductivity is σ a ¼ σ 1 CSðz1 Þ þ

N −1 X

σ i ½CSðzi Þ−CSðzi−1 Þ þ σ N ½1−CSðzN−1 Þ;

ðA2Þ

i¼2

where σi and zi(i = 1,2,…,N) are layers' conductivity and normalized depth to the interfaces, respectively. Appendix B Koefoed (1979) developed algorithm to forward model resistivity response of 1-D earth models with a four-electrode system: R ¼ Q ðjAα jÞ−Q ðjBα jÞ−Q ðjAβjÞ þ Q ðjBβjÞ;

ðB1Þ

with A, B, α, and β indicating electrode locations and, ∞

Q ðrÞ ¼

ρ1 ∫ K 1 ðλÞJ 0 ðλr Þdλ; 2π

ðB2Þ

0

where J0 is Bessel function of order zero and λ is the variable of integration and

K n ðλÞ ¼

K nþ1 þ σ nþ1 σn

σ nþ1 σi

tanhðλhn Þ

þ K nþ1 tanhðλhn Þ

; K N ¼ 1;

ðB3Þ

where σn and hn are resistivity and thickness of layer n. IngemanNielsen and Baumgartner (2006) indicate that the computational efficiency of Eq. (B1) may be increase by using ( ) ∞ ρ1 1 Q ðrÞ ¼ ; ∫ ½K ðλÞ−1 J 0 ðλr Þdλ þ 2π 0 1 r

ðB4Þ

which is obtained by subtracting the homogeneous half-space response from the kernel and adding it outside the Hankel transform.

53

References Agostinetti, N.P., Malinverno, A., 2010. Receiver function inversion by trans-dimensional Monte Carlo sampling. Geophysical Journal International 181, 858–872. Ahmed, N.A., Gokhale, D.V., 1989. Entropy expressions and their estimators for multivariate distributions. IEEE Transaction on Information Theory 35 (3), 688–692. Atekwana, E.A., Atekwana, E.A., 2010. Geophysical signatures of microbial activity at hydrocarbon contaminated sites: a review. Surveys in Geophysics 31, 247–283. Auken, E., Christiansen, A.V., Jacobsen, B.H., Foged, N., Sørensen, K.I., 2005. Piecewise 1D laterally constrained inversion of resistivity data. Geophysical Prospecting 53, 497–506. Bedard, M., 2008. Optimal acceptance rates for Metropolis algorithms: moving beyond 0.234. Stochastic Processes and Their Applications 118 (12), 2198–2222. http:// dx.doi.org/10.1016/j.spa.2007.12.005. Binley, A., Kemna, A., 2005. Electrical methods. In: Rubin, Y., Hubbard, S. (Eds.), Hydrogeophysics. Springer, pp. 129–156. Bodin, T., Sambridge, M., 2009. Seismic tomography with the reversible jump algorithm. Geophysical Journal International 178, 1411–1436. Bodin, T., Sambridge, M., Rawlinson, N., Arroucau, P., 2012. Transdimensional tomography with unknown data noise. Geophysical Journal International. http://dx.doi.org/ 10.1111/j.1365-246X.2012.05414.x. Bogaert, P., Fasbender, D., 2007. Bayesian data fusion in a spatial prediction context: a general formulation. Stochastic Environmental Research Risk Assessment 21 (6), 695–709. Brodie, R., Sambridge, M., 2006. A holistic approach to inversion of frequency-domain airborne EM data. Geophysics 71 (6), G301–G312. Brosten, T.R., Day-Lewis, F.D., Schultz, G.M., Curtis, G.P., Lane Jr., J.W., 2011. Inversion of multi-frequency electromagnetic induction data for 3D characterization of hydraulic conductivity. Journal of Applied Geophysics 73 (4), 323–335. Cameron, D.R., De Jong, E., Read, D.W.L., Oosterveld, M., 1981. Mapping salinity using resistivity and electromagnetic techniques. Canadian Journal of Soil Sciences 61, 67–78. Chen, J., Hoversten, G.M., Vasco, D., Rubin, Y., Hou, Z., 2007. A Bayesian model for gas saturation estimation using marine seismic AVA and CSEM data. Geophysics 72, WA85–WA95. Christakos, G., 1990. A Bayesian/maximum-entropy view to the spatial estimation problem. Mathematical Geology 22 (7), 763–777. Christakos, G., 2000. Modern Spatiotemporal Geostatistics. Oxford University Press, New York, NY (304 pp.). Clemen, R.T., Winkler, R.L., 1993. Aggregating point estimates: a flexible modeling approach. Management Sciences 39, 501–515. Commer, M., Newman, G.A., 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion. Geophysical Journal International 178 (3), 1305–1316. Cover, T.M., Thomas, J.A., 2006. Elements of Information Theory. Wiley. Cowles, M.K., Carlin, B.P., 1996. Markov chain Monte Carlo convergence diagnostics: a comparative review. Journal of the American Statistical Association 91 (434), 883–904. Cressie, N.A.C., 1985. Fitting Variogram Models by Weighted Least Squares. Mathematical Geology 17, 563–586. Crook, N., Binley, A., Knight, R., Robinson, D.A., Zarnetske, J., Haggerty, R., 2008. Electrical Resistivity Imaging of the Architecture of Sub-Stream Sediments. Water Resources Research 44, W00D13. http://dx.doi.org/10.1029/2008WR006968. Curtis, A., 2004a. Theory of model-based geophysical survey and experimental design Part A Linear Problems. The Leading Edge 23 (10), 997–1004. Curtis, A., 2004b. Theory of model-based geophysical survey and experimental design Part B Nonlinear Problems. The Leading Edge 23 (10), 1112–1117. Doetsch, J.A., Linde, N., Coscia, I., Greenhalgh, S., Green, A.G., 2010a. Zonation for 3D aquifer characterization based on joint inversion of multi-method crosshole geophysical data. Geophysics 75 (6), G53–G64. Doetsch, J., Linde, N., Binley, A., 2010b. Structural joint inversion of time-lapse crosshole ERT and GPR traveltime data. Geophysical Research Letters 37, L24404 (ISSN 0094–8276). Douaik, A., Van Meirvenne, M., To´th, T., Serre, M.L., 2004. Space–time mapping of soil salinity using probabilistic Bayesian maximum entropy. Stochastic Environmental Research Risk Assessment 18, 219–227. Dubreuil-Boisclair, C., Gloaguen, E., Marcotte, D., Giroux, B., 2011. Heterogeneous aquifer characterization from GPR tomography and borehole hydrogeophysical data using nonlinear Bayesian simulations. Geophysics 76 (4), 1–13. Eni, 2008. Attiva di monitoraggio nell'area oggetto degli intervent di bonifica a seguito dell’incidente al pozzo, TR 24, ENI S.p.A. Fasbendar, D., Brasseur, O., Bogaert, P., 2009. Bayesian data fusion for space–time prediction of air pollutants: The case of NO2 in Belgium. Atmospheric Environment 43, 4632–4645. Gallardo, L., Meju, M., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross gradient constraints. Journal of Geophysical Research 109, 1–11. Gelman, A., Rubin, D., 1992. Inference from iterative simulation using multiple sequences. Statistical Science 7, 457–472. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian Data Analysis. C. Hall, Chapman and Hall (526 pp.). Gelman, A., Roberts, G.O., Gilks, W.R., 1996. Efficient Metropolis jumping rules. In: Bernardo, J.M., Berger, J.O., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics, 5. Oxford University Press, pp. 599–607. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1996. Markov Chain Monte Carlo in Practice. Chapman and Hall/CRC Press, London.

54

A. JafarGandomi, A. Binley / Journal of Applied Geophysics 96 (2013) 38–54

Green, P.J., 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82, 711–732. Green, P., Mira, A., 2001. Delayed rejection in reversible jump Metropolis-Hastings. Biometrika 88 (4), 1035–1053. Guest, T., Curtis, A., 2009. Iteratively constructive sequential design of experiments and surveys with nonlinear parameter-data relationships. Journal of Geophysical Research 114, B04307. http://dx.doi.org/10.1029/2008JB005948. Guest, T., Curtis, A., 2010. Optimal trace selection for amplitude versus angle (AVA) processing of shale-sand reservoirs. Geophysics 75 (4), C37–C47. Gyulai, A´., Ormos, T., 1999. A new procedure for the interpretation of VES data: 1.5-D simultaneous inversion method. Journal of Applied Geophysics 64, 1–17. Hastings, W., 1970. Monte Carlo simulation methods using Markov chains and their applications. Biometrika 57, 97–109. Hezarjaribi, A., Sourell, H., 2007. Feasibility study of monitoring the total available water content using non-invasive electromagnetic induction-based and electrodebased soil electrical conductivity measurements. Irrigation and Drainage 56, 53–65. Hoversten, G.M., Cassassuce, F., Gasperikova, E., Newman, G.A., Chen, J., Rubin, Y., Hou, Z., Vasco, D., 2006. Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data. Geophysics 71 (3), C1–C13. Hu, W., Abubakar, A., Habashy, T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints. Geophysics 74 (6), R99–R109. Huisman, J.A., Hubbard, S.S., Redman, J.D., Annan, P., 2003. Measuring soil water content with ground penetrating radar. Vadose Zone Journal 2 (4), 476–491. Ingeman-Nielsen, T., Baumgartner, F., 2006. CR1Dmod: a Matlab program to model 1D complex resistivity effects in electrical and electromagnetic surveys. Computers and Geosciences 32, 1411–1419. JafarGandomi, A., Curtis, A., 2012. Assessing the monitorability of CO2 saturation in subsurface aquifers. International Journal of Greenhouse Gas Control 7, 244–260. http://dx.doi.org/10.1016/j.ijggc.2011.10.015. Jouini, M.N., Clemen, R.T., 1996. Copula models for aggregating expert opinions. Operations Research 44, 444–457. Kaufman, A.A., Keller, G.V., 1983. Frequency and transient soundings. Methods in Geochemistry and Geophysics, 16. Elsevier, NY. Koefoed, O., 1979. Geosounding Principles, 1. Methods in Geochemistry and Geophysics, 14A. Elsevier, Amsterdam (276 pp.). Kosinski, W.K., Kelly, W.E., 1981. Geoelectric soundings for predicting aquifer properties. Ground Water 19 (2), 163–171. Krause, A., Singh, A., Guestrin, C., 2008. Near-optimal sensor placements in Gaussian processes: theory, efficient algorithm and empirical studies. Journal of Machine Learning Research 9, 235–284. Krykacz-Haussmann, B., 2001. Epistemic sensitivity analysis based on the concept of entropy. In: Prado, P., Bolado, R. (Eds.), Proceedings of SAMO2001. CIEMAT, Madrid, pp. 31–35. Kvamme, K.L., 2006. Integrating multidimensional geophysical data. Archaeological Prospection 13 (1), 57–72. Kvamme, K.L., 2007. Integrating multiple geophysical datasets. In: Wiseman, J., El-Baz, F. (Eds.), Remote sensing in archaeology. Springer, New York, pp. 345–374. Linde, N., Doetsch, J.A., 2010. Joint Inversion of Crosshole GPR and Seismic Traveltime Data. In: Miller, R.D., Bradford, J.H., Holliger, K. (Eds.), Advances in Near-surface Seismology and Ground-penetrating Radar. Society of Exploration Geophysicists, pp. 1–18. http://dx.doi.org/10.1190/1.9781560802259.ch1 (Ch. 1). Lindley, D., 1956. On a measure of the information provided by an experiment. Annals of Mathematical Statistics 27, 986–1005. Looms, M.C., Binley, A., Jensen, K.H., Nielsen, L., Hansen, T.M., 2008. Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data. Vadose Zone Journal 7 (1), 238–248. Malinverno, A., 2002. Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem. Geophysical Journal International 151, 675–688. Maurer, H., Curtis, A., Boerner, D., 2010. Recent advances in optimized geophysical survey design. Geophysics 75 (5) (75A177-75A194). McClymont, A.F., Roy, J.W., Hayashi, M., Bentley, L.R., Maurer, H., Langston, G., 2011. Investigating groundwater flow paths within proglacial moraine using multiple geophysical methods. Journal of Hydrology 399, 57–69.

McNeill, J.D., 1980. Electromagnetic terrain conductivity measurement at low induction numbers. Tech. Note TN-6.Geonics Ltd., Mississauga, ON. Meles, G.A., Greenhalgh, S.A., Green, A.G., Maurer, H., Van der Kruk, J., 2012. GPR full waveform sensitivity and resolution analysis using an FDTD adjoint method. IEEE Transactions on Geosciences and Remote Sensing 50 (5), 1881–1896. Metropolis, N., Rosenbluth, M.N., Rosenbluth, A.W., Teller, A.H., Teller, E., 1953. Equation of state calculations by fast computing machines. Journal of Chemical Physics 21, 1087–1092. Minsley, B.J., 2011. A trans-dimensional Bayesian Markov chain Monte Carlo algorithm for model assessment using frequency-domain electromagnetic data. Geophysical Journal International 187, 252–272. Oldenburg, D.W., Li, Y., 1994. Inversion of induced polarization data. Geophysics 59, 1327–1341. Pride, S., 1994. Governing equations for the coupled electromagnetics and acoustic of porous media. Physical Review B50 (21), 15678–15696. Roberts, G.O., Gelman, A., Gilks, W.R., 1997. Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability 7, 110–120. Rovetta, D., Lovatini, A., Watts, M.D., 2008. Probabilistic joint inversion of TD-CSEM, MT and DC data for hydrocarbon exploration. SEG Expanded Abstracts, 27, p. 599. http://dx.doi.org/10.1190/1.3063723. Sambridge, M., Gallagher, K., Jackson, A., Rickwood, P., 2006. Transdimensional inverse problems, model comparison and the evidence. Geophysical Journal International 167, 528–542. Santos, F.A.M., 2004. 1-D laterally constrained inversion of EM34 profiling data. Journal of Applied Geophysics 56, 123–134. Sebastiani, P., Wynn, H.P., 2000. Maximum entropy sampling and optimal Bayesian experimental design. Journal of the Royal Statistical Society 62, 145–157. Serre, M.L., 1999. Environmental Spatiotemporalmapping and Ground Water Flow Modelling Using the BME and ST Methods. University of North Carolina, Chapel Hill, NC. Serre, M.L., Christakos, G., 1999. Modern geostatistics: computational BME in the light of uncertain physical knowledge—the Equus Beds study. Stochastic Environmental Research and Risk 13 (1), 1–26. Serre, M.L., Bogaert, P., Christakos, G., 1998. Computational investigations of Bayesian maximum entropy spatiotemporal mapping. In: Buccianti, A., Nardi, G., Potenza, R. (Eds.), 4th Annual Conference, vol. 1. De Frede Editore, Naples, Italy, pp. 117–122. Shannon, C.E., 1948. A mathematical theory of communication. Bell Labs Technical Journal 27, 623–656. Sivia, D., 1996. Data Analysis: A Bayesian Tutorial. Oxford University Press, USA. Smith, J.T., Booker, J.R., 1991. Rapid inversion of two and threedimensional magnetotelluric data. Journal of Geophysical Research 96, 3905–3922. Tarantola, A., 2005. Inverse Problem Theory and Model Parameter Estimation. Society for Applied and Industrial Mathematics, Philadelphia. Tierney, L., Mira, A., 1999. Some adaptive Monte Carlo methods for Bayesian inference. Statistics in Medicine 18, 2507–2515. Titov, K., Revil, A., Konosavsky, P., Straface, S., Troisi, S., 2005. Numerical modelling of self-potential signals associated with a pumping test experiment. Geophysical Journal International 162, 641–650. Triantafilis, J., Monteiro Santos, F.A., 2011. Hydrostratigraphic analysis of the Darling River valley (Australia) using electromagnetic induction data and a spatially constrained algorithm for quasi-three-dimensional electrical conductivity imaging. Hydrogeology Journal 19, 1053–1063. van den Berg, J., Curtis, A., Trampert, J., 2003. Optimal nonlinear Bayesian experimental design: an application to amplitude versus offset experiments. Geophysical Journal International 155, 411–421. Viezzoli, A., Christiansen, A.V., Auken, E., Sørensen, K., 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion. Geophysics 73, F105–F113. Vrugt, J.A., ter Braak, C.J.F., Diks, C.G.H., Robinson, B.A., Hyman, J.M., Higdon, D., 2009. Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. International Journal of Nonlinear Sciences and Numerical Simulation 10 (3), 273–290. Zidek, J.V., Sun, W., Le, N.D., 2000. Designing and integrating composite networks for monitoring multivariate Gaussian pollution fields. Applied Statistics 49 (1), 63–79.