Usually expressed as a percentage, lag tolerance corresponds to how much the distance between sample pairs can differ from the exact lag distance and still be included in the lag calculations. Tolerance of 0.5 or 50% is usually chosen for the calculation of the lag tolerance.

Figure 76. Example of an isotropic variogram; lag = 20 m; tolerance lag = 0.5, using Isatis.

Angle

Environmental data often has different levels of spatial correlation in different directions, which is referred to as anistropy. It is essential to account for anisotropy in variograms. Generally, the variogram map is used as a visualization tool for identifying anisotropic behaviors (see Figure 77). Anisotropy can be quantified by examining variograms constructed from data pairs along different directions. For example, it is possible to build a variogram by dividing the plane into four equal parts of 45°. Then, four variograms are constructed: one for the direction where the maximum spatial correlation is observed (reference direction) and three other perpendicular directions in order to define the lower levels of correlation in other directions. It is also possible to calculate the short-range variogram by taking shorter lags in order to specify the behavior at the origin.

Figure 77. Example of an anisotropic variogram taking into account the anisotropic direction (N110°) as reference, using Isatis.

Figure 78. Variogram cloud (left) and experimental variogram (isotropic, log value = 20 m) (right), using Isatis.

For example, if contaminant concentration varies rapidly, the initial gradient of the empirical variogram will be steep. In this case, an exponential or linear variogram may better represent this behavior. If a noncontinuous behavior is expected (inorganic contamination, for instance) or if analytical or sampling error is suspected, the nugget effect may represent such uncertainties.

Theoretical variogram models commonly used include exponential, spherical, Gaussian, Matérn, and linear. With the exception of the linear model, these theoretical models differ primarily in their behavior close to the origin. An illustration of the common theoretical variogram models is shown in Figure 79.

Figure 79. Theoretical variogram models.

When fitting a variogram, whether to use a nugget effect, and the most suitable range or sill value must be determined. It is best to produce several variogram models that could match with the contaminant concentration behavior. Then, cross validation can be used to compare the accuracy of each model in order to identify the most appropriate.

In the cross-validation process, several errors can be calculated (such as mean standardized error or mean error). A robust model shows a standard error of the variance close to 1, a strong correlation between the true values and the estimated values and, finally, a standard error close to 0. Through the analysis of these parameters, the most suitable model (the most accurate) can be chosen. Figure 80 illustrates a cross-validation example and error analysis.

Figure 80. Cross-validation example and error analysis using Isatis.

Kriging is known as the best linear unbiased estimator because:

• The search neighborhood weights are chosen so that the estimator is unbiased (difference between the predicted value and measured realization equals zero).
• Kriging solves a series of linear equations (kriging system) by minimizing the kriging variance.

Different kriging methods are unique in their assumptions, constraints, limitations, and purpose. Popular kriging methods include ordinary kriging, indicator kriging, simple kriging (see conditional simulation), factorial kriging, and universal kriging. When more than one variable contributes to the kriging estimate in each of these methods, it is termed “co-kriging.”

Kriging can also be used to upscale (estimate over a given area or volume) or downscale (estimate at a specific location). Block kriging is a form of upscaling spatial information across a sampling domain. Point, or punctual, kriging is a means for downscaling spatial information across a sampling domain. The use of point or block kriging can be guided by the relationship of the sample design (support and spacing) to the interpolation grid cell size.

Often, many different types of data are collected from a contaminated site to directly and indirectly assess factors controlling the fate and distribution of contaminants. In soil systems, for example, it is common to integrate environmental geophysics, LiDAR, and soil sample information to characterize a site in space over time. Measured parameters that exhibit shared structured spatial variation (spatial autocorrelation or shared space-dependent variation) can be used together for interpolation, which is the process of co-kriging.

The search neighborhood is optimized through the cross variogram for co-kriging. Co-kriging is popular for integrating a high-resolution data set with a more sparsely sampled data set to generate spatial estimates with greater coverage. Co-kriging can be applied as block co-kriging or point co-kriging.

Collected data often has a skewed distribution, which can introduce significant bias in the geospatial estimates, for example, when using ordinary kriging. In this case, the mean of the distribution is not representative of the data set. Ordinary kriging, also called the mean interpolator, performs the estimates of concentrations by using local means and the estimated values will be closely calculated around this mean. When an asymmetric distribution is observed, the kriged data histogram will not be consistent with the histogram of the original data.

Performing logarithmic or Gaussian (normalizing) transformations of the raw distributions may remediate this effect. If they are being used, data transformation should be completed before proceeding with kriging methods. Moreover, Gaussian (normal) distributions are needed to perform the conditional simulation method. This type of distribution may minimize the influence of high values compared to the low ones, and reveal the best data distribution in order to obtain a more structured experimental variogram

Below are general descriptions of point kriging and block kriging. These kriging approaches can be applied to all kriging methods discussed in this section.

Point kriging is one of the most common forms of kriging. Point kriging generates a point-estimated value of a sampled property at an unsampled location, using optimally weighted, known (sampled) ambient values. Point kriging is useful for estimating or simulating a measured property at finer spatial resolutions than was used in sampling (downscaling). Point kriging estimates are sensitive to local discontinuities and large nugget effects.

Block kriging estimates a weighted average across a particular domain, in this case a grid cell or block. A predefined, regular-spaced interpolation grid is necessary to perform block kriging, which is ideal for studying regional patterns of variation, but not for local-scale variation. Block kriging inherently generates smoother maps than point kriging, helping to negate the effects of discontinuities in the data as well as other undesirable artifacts. The block kriging variance is lower in comparison to point kriging because the local scale variation, or within-block variation, is removed—an advantage for data sets with a large nugget.

Simple kriging is used for spatial interpolation when the mean and spatial correlation model are constant and known. Simple kriging is applied in conditional simulation.

This method is subject to the following constraints:

• Data are second-order stationary if covariance function model is used or intrinsically stationary if variogram model is used.
• Spatial correlation is present among the data.
• No duplicate sites are present.

The spatial correlation model for kriging can be based on either the covariance function or variogram. In general, it is easier to estimate a variogram model than a covariance function because no estimate of the mean is required.

The pure nugget effect indicates the absence of spatial correlation. In this case, a deterministic geospatial method should be considered.

This method assumes the mean and variance are constant and known (second-order stationarity).

This method is very restrictive due to the requirement that the mean is known.

In addition to interpolated values, kriging can provide an estimate of the uncertainty in the interpolated values, which is known as the kriging variance or standard error. The quality of the kriging model fit to the data should be evaluated using cross validation or validation.

In the context of optimization, see how to use the results of the geospatial methods to address specific optimization questions.

Ordinary kriging is used to generate estimates of a sampled variable in unsampled locations.

This method is subject to the following constraints:

• Intrinsic stationarity—the variogram has a sill.
• No duplicate sites are used.
• The nugget is small relative to the sill.

The spatial correlation model for kriging can be based on either the covariance function or variogram. In general, it is easier to estimate a variogram model than a covariance function because no estimate of the mean is required.

A pure nugget effect observed in the variogram indicates the absence of spatial correlation (strict stationarity). A variogram with a large nugget relative to the sill indicates variability in the data that cannot be accurately predicted using kriging. The nugget effect and shape of the variogram near the origin have the most influence on the kriging predictions. Since observations located beyond the range indicated by the variogram are spatially uncorrelated, the range can be used to select the kriging search radius.

• The data are locally stationary within the kriging search neighborhood.
• Data are locally second-order stationary if covariance function is used or intrinsically stationary if variogram is used.
• The data are spatially dependent.
• The expected value of a sample is the mean of the population.
• The difference in values between two samples is due to the distance between the two locations.

This kriging method is flexible enough to give good results for most data sets, unless there are strong spatial trends.

In addition to interpolated values, kriging can provide an estimate of the uncertainty in the interpolated values, which is known as the kriging variance or standard error. The quality of the kriging model fit to the data should be evaluated using cross validation or validation.

In the context of optimization, see how to use the results of the geospatial methods to address specific optimization questions.

Universal kriging and kriging with an external trend are used for spatial interpolation when there is a trend in the data. In applications for remediation optimization, this approach may include distribution of contaminant concentrations from a single point source.

It is necessary to develop a model of the semivariogram or the covariance function. Several software packages can be used to perform the necessary operations, with varying requirements for user expertise.

• Data have a spatial trend that is modeled using a simple function of the location coordinates.
• Detrended data are locally second-order stationary if covariance function is used or intrinsically stationary if variogram is used.

Universal kriging provides the ability to interpolate spatially related values when there is a trend in the mean, a property not realized with ordinary kriging. However, the use of universal kriging is limited to scenarios in which the variation exhibited by the attribute of interest can be considered to be systematic. As with most interpolation techniques, the uncertainty associated with the results of the method increases dramatically with the complexity of the model used, and should be taken into account.

In addition to interpolated values, kriging can provide an estimate of the uncertainty in the interpolated values known as the kriging variance or standard error. The quality of the kriging model fit to the data should be evaluated using cross validation or validation.

In the context of optimization, see how to use the results of the geospatial methods to address specific optimization questions.

This method is typically used to map the probabilities of exceeding the specified threshold.

The concentrations must be transformed into a binary variable (0,1). Several variograms must be developed and fitted if several thresholds are set.

The verifications that must be performed in order to implement this method are:

• the coherence between the variograms constructed for each indicator
• the cross validation for each variogram model and the neighborhood (taken into account for the interpolation)
• the variance map

This method assumes stationarity of each indicator.

The advantages of indicator kriging method include:

• The method does not require any assumption about data distribution.
• The method can be used for categorical data (type of soil, index of contamination).

The weaknesses of this method are:

• The same value (1) is assigned for all the sample points exceeding the specified threshold. The degree of contamination is not taken into account. This means a result several orders of magnitude above the threshold has the same effect on nearby exceedance probability as a result that is less than 1 percent above the threshold.
• A variogram must be constructed and fitted for each indicator corresponding to each threshold.
• The use of a nugget effect is questionable for an application where it is assumed that comparisons to a threshold value are exact (Krivoruchko and Bivand 2009).
• Negative probabilities or probabilities greater than 1 may appear because of the weights applied by the kriging or if the variograms for each threshold are not compatible.

A major limitation to kriging techniques and their associated estimates of exceedance probability is that these methods minimize the variance to generate spatial estimates. This result inherently underestimates the variability of the system being characterized and thereby biases the probability outcomes. An alternative to kriging that overcomes this limitation is conditional simulation.

The results obtained by performing indicator kriging are expressed in terms of probability. Thus, mapping of the results does not indicate the concentrations of the contaminants, but rather the risk of exceeding a threshold such as the remediation goal.

In the context of optimization, see how to use the results of the geospatial methods to address specific optimization questions.

Probability kriging and disjunctive kriging are two alternatives that can overcome some of the limitations of indicator kriging. Probability kriging uses autocorrelation of both the original data and the indicator binary variable, and cross-correlation between them. Probability kriging uses co-kriging to preserve some of the information lost in indicator kriging (ESRI 2013).

Disjunctive kriging is defined as the simple kriging of data transformed to a standard normal distribution using Hermite polynomials. A probability distribution and associated map of probabilities can then be generated by comparing estimated values with the normal distribution (Webster and Oliver 2001). Disjunctive kriging assumes that the data in question follow a bivariate normal distribution. This assumption means that linear combinations of paired values are normally distributed with correlation coefficients depending solely on separation distance rather than absolute position (Krivoruchko 2011). Probability kriging is available in ArcGIS Geostatistical Analyst, and disjunctive kriging is available in ArcGIS Geostatistical Analyst, Isatis, and R.

Regression kriging is a hybrid of regression and kriging methods. A regression to predict the quantity of interest based on auxiliary variables (such as land surface elevation or soil type) is performed first, and then simple kriging is performed on the regression residuals. This method is mathematically equivalent to kriging with external drift (trend), where the auxiliary variable relationship and kriging model is solved simultaneously to calculate the kriging weights. The advantage of regression kriging is that it can extend the method to a broader range of regression techniques and allow separate interpretation of the two interpolated components.

Factorial Kriging Analysis (FKA) involves multivariate variogram modeling, principal component analysis, and co-kriging. Factorial kriging allows investigators to simultaneously study spatial relationships between multiple variables according to their shared spatial ranges of autocorrelation. Three primary outputs are generated from this technique:

• regionalized correlation coefficients
• regionalized factors (principal components)
• individual sets of concentration maps (univariate) and principal component maps (multivariate)

Spatial autocorrelation implies that a fraction of the total measured variation in the data are space dependent, or structured. Although several analytical tools can model spatial autocorrelation, the variogram is the fundamental component of factorial kriging.

Integrating densely sampled secondary measurements can lead to more consistent descriptions of sparsely sampled direct measurements (Castrignanò et al. 2012). When the secondary variable is sampled much more densely than the primary variable, however, it can cause instability in solving the kriging system because the correlation between close secondary data is greater than that of sparse primary data (Goovaerts 1997). Collocated co-kriging is an efficient method to integrate exhaustive secondary measurements with sparse primary measurements to estimate the primary variable at target locations (Morari, Castrignanó, and Pagliarin 2009). Collocated co-kriging is similar to ordinary co-kriging, except that the neighborhood search specifically uses the secondary variable information collocated with the measured primary variable and the target location to be estimated. In depth discussion of collocated co-kriging is beyond the scope of this document; suggested readings on this topic include Goovaerts 1997; Castrignanò et al. 2012; and Morari, Castrignanó, and Pagliarin 2009.

The results obtained by performing conditional simulation have several applications:

• mapping the contamination sources according to a remediation goal
• mapping the uncertainties associated with the estimates in order to determine the financial risk of the project
• defining new sampling programs according to the mapped uncertainties
• defining the strategy for cleanup

This method requires a Gaussian transformation of the original data set. Several techniques allow the transformation of the raw data to a Gaussian distribution (for example, Hermite polynomials, Gaussian anamorphosis). In order to produce a conditional simulation, a new variogram must be developed with the Gaussian distribution and then fitted. Variogram development and simulations must be performed on Gaussian transformed data and results must be back-transformed.

This method assumes normal distribution and second order stationarity if simple kriging is used.

The advantages of the conditional simulation method are:

• Each simulation respects the values observed in the sample points. Mean and variance are preserved in each realization.
• Performing a large number of simulations better reproduces the real diversity of possible cases, when random variables are modeled.
• It is possible to quantify and locate the areas of uncertainty and to support management of the financial risks of a project.

The weaknesses are:

• Some distributions do not allow a perfect Gaussian transformation of the data set.
• The time for completion can be longer than other methods.
• Because of the large data sets generated, computation time may be significant.

The results of conditional simulation are expressed by a number of equally probable maps. Each realization obtained by using the conditional simulation method is equally valid (Figure 81).

Figure 81. Conditional simulation results.

Which is the most appropriate model to represent the contamination? A simulation postprocessing step compiles the set of models produced by this method. The objective of this postprocessing is to answer the following questions:

1. A) What are the contaminated areas? This question is used to locally estimate the risk of exceeding a remediation threshold.

When a remediation threshold is applied, a first postprocessing can be performed in order to estimate the local level of contamination within each block. For example, a block x is simulated a hundred times. If this block x exceeds the remediation threshold 20 times, then the local probability of exceeding the threshold for this block x is 20%. This process is repeated for each block of the grid.

The final maps are expressed in terms of probability which defines the contamination risk within each block according to the remediation threshold. Figure 82 shows the results.

Figure 82. Probability map.

1. B) What is the contaminated volume? This question is used to estimate the overall level of contamination in the study area.

Similarly, the postprocessing of conditional simulation allows you to estimate the overall degree of contamination on the site. If a hundred simulations are performed, the contaminated volume can be calculated for each realization, by calculating the number of blocks showing concentrations above the remediation threshold. Thus, if a hundred simulations are performed, a hundred volumes of contamination will be obtained. This technique allows you to define the overall degree of contamination on the site and to define a margin of error by obtaining the optimistic and pessimistic volume (percentiles 5 and 95) and the most probable volume (median of the distribution). Figure 83 shows the results of this technique.

Figure 83. Using conditional simulation results to determine probable volume of contamination.

An accurate estimate minimizes the differences between the optimistic and the pessimistic volumes.

Furthermore, other maps can be produced to better interpret the results:

• mean of the realizations: a map similar to the results obtained by kriging
• smallest and largest realizations: shows the optimistic and the pessimistic mappings obtained by performing the conditional simulation
• standard deviation of realizations: quantifies the error of the estimation within each block and gives an estimate of the calculation accuracy

The verifications that must be performed in order to implement this method are:

• Gaussian transformation of the data set
• cross-validation for the variogram model and the neighborhood taken into account for the interpolation
• quantification of the uncertainties and the differences between the optimistic and pessimistic volumes in order to determine the measure of certainty of the estimates.

In the context of optimization, see how to use the results of the geospatial methods to address specific optimization questions.

Co-Simulation

Co-simulation is an advanced method that is used for estimating values for a primary variable when there are secondary data available (known as proxy data). It is similar to conditional simulation. In co-simulation multiple values are generated for the primary variable at unsampled locations using both the primary and secondary sampled data. The values are generated using a probability simulation technique. Like co-kriging, co-simulation is used when the primary variable is under-sampled and the secondary data are spatially correlated with the primary variable. A cross variogram is used to check the spatial correlation; see co-simulation for estimating concentrations based on proxy data.

• “PAH contamination in sediments: Uncertainty analysis,” case study that uses conditional simulation (Section 6.2)
• “The turning bands method for simulation of random fields using line generation by a spectral method” (Mantoglou and Wilson 1982)
• Geostatistics: Modeling Spatial Uncertainty, Second Edition (Chilès and Delfiner 2012)
• Geostatistical Simulation, Models and Algorithms (Lantuéjoul 2002)
• Plurigaussian Simulations in Geosciences (Armstrong et al. 2011)
• Introduction to Disjunctive Kriging and Nonlinear Geostatistics (Rivoirard 1994)