Geostatistical analysis and conditional simulation for estimating the spatial variability of hydraulic conductivity in the Choushui River alluvial fan,Taiwan

This work evaluated the spatial variability and distribution of heterogeneous hydraulic conductivity (K) in the Choushui River alluvial fan in Taiwan, using ordinary kriging (OK) and mean and individual sequential Gaussian simulations (SGS). A baseline flow model constructed by upscaling parameters was inversely calibrated to determine the pumping and recharge rates. Simulated heads using different K realizations were then compared with historically measured heads. A global/local simulated error between simulated and measured heads was analysed to assess the different spatial variabilities of various estimated K distributions. The results of a MODFLOW simulation indicate that the OK realization had the smallest sum of absolute mean simulation errors (SAMSE) and the SGS realizations preserved the spatial variability of the measured K fields. Moreover, the SAMSE increases as the spatial variability of the K field increases. The OK realization yields small local simulation errors in the measured K field of moderate magnitude, whereas the SGS realizations have small local simulation errors in the measured K fields, with high and low values. The OK realization of K can be applied to perform a deterministic inverse calibration. The mean SGS method is suggested for constructing a K field when the application focuses on extreme values of estimated parameters and small calibration errors, such as in a simulation of contaminant transport in heterogeneous aquifers. The individual SGS realization is useful in stochastically assessing the spatial uncertainty of highly heterogeneous aquifers. Copyright © 2004 John Wiley & Sons, Ltd.     

A comparative study of Monte Carlo simple genetic algorithm and noisy genetic algorithm for cost-effective sampling network design under uncertainty

This study evaluates and compares two methodologies, Monte Carlo simple genetic algorithm (MCSGA) and noisy genetic algorithm (NGA), for cost-effective sampling network design in the presence of uncertainties in the hydraulic conductivity (K) field. Both methodologies couple a genetic algorithm (GA) with a numerical flow and transport simulator and a global plume estimator to identify the optimal sampling network for contaminant plume monitoring. The MCSGA approach yields one optimal design each for a large number of realizations generated to represent the uncertain K-field. A composite design is developed on the basis of those potential monitoring wells that are most frequently selected by the individual designs for different K-field realizations. The NGA approach relies on a much smaller sample of K-field realizations and incorporates the average of objective functions associated with all K-field realizations directly into the GA operators, leading to a single optimal design. The efficacy of the MCSGA-based composite design and the NGA-based optimal design is assessed by applying them to 1000 realizations of the K-field and evaluating the relative errors of global mass and higher moments between the plume interpolated from a sampling network and that output by the transport model without any interpolation. For the synthetic application examined in this study, the optimal sampling network obtained using NGA achieves a potential cost savings of 45% while keeping the global mass and higher moment estimation errors comparable to those errors obtained using MCSGA. The results of this study indicate that NGA can be used as a useful surrogate of MCSGA for cost-effective sampling network design under uncertainty. Compared with MCSGA, NGA reduces the optimization runtime by a factor of 6.5.     

Estimating threshold-exceeding probability maps of environmental variables with Markov chain random fields

Estimating and mapping spatial uncertainty of environmental variables is crucial for environmental evaluation and decision making. For a continuous spatial variable, estimation of spatial uncertainty may be conducted in the form of estimating the probability of (not) exceeding a threshold value. In this paper, we introduced a Markov chain geostatistical approach for estimating threshold-exceeding probabilities. The differences of this approach compared to the conventional indicator approach lie with its nonlinear estimators—Markov chain random field models and its incorporation of interclass dependencies through transiograms. We estimated threshold-exceeding probability maps of clay layer thickness through simulation (i.e., using a number of realizations simulated by Markov chain sequential simulation) and interpolation (i.e., direct conditional probability estimation using only the indicator values of sample data), respectively. To evaluate the approach, we also estimated those probability maps using sequential indicator simulation and indicator kriging interpolation. Our results show that (i) the Markov chain approach provides an effective alternative for spatial uncertainty assessment of environmental spatial variables and the probability maps from this approach are more reasonable than those from conventional indicator geostatistics, and (ii) the probability maps estimated through sequential simulation are more realistic than those through interpolation because the latter display some uneven transitions caused by spatial structures of the sample data.     

Plurigaussian modeling of geological domains based on the truncation of non-stationary Gaussian random fields

The plurigaussian model is used in mining engineering, oil reservoir characterization, hydrology and environmental sciences to simulate the layout of geological domains in the subsurface, while reproducing their spatial continuity and dependence relationships. However, this model is well-established only in the stationary case, when the spatial distribution of the domains is homogeneous in space, and suffers from theoretical and practical impediments in the non-stationary case. To overcome these limitations, this paper proposes extending the model to the truncation of intrinsic random fields of order k with Gaussian generalized increments, which allows reproducing spatial trends in the distribution of the geological domains. Methodological tools and algorithms are presented to infer the model parameters and to construct realizations of the geological domains conditioned to existing data. The proposal is illustrated with the simulation of rock type domains in an ore deposit in order to demonstrate its applicability. Despite the limited number of conditioning data, the results show a remarkable agreement between the simulated domains and the lithological model interpreted by geologists, while the conventional stationary plurigaussian model turns out to be unsuccessful.     

Properties and limitations of sequential indicator simulation

The sequential indicator algorithm is a widespread geostatistical simulation technique that relies on indicator (co)kriging and is applicable to a wide range of datasets. However, such algorithm comes up against several limitations that are often misunderstood. This work aims at highlighting these limitations, by examining what are the conditions for the realizations to reproduce the input parameters (indicator means and correlograms) and what happens with the other parameters (other two-point or multiple-point statistics). Several types of random functions are contemplated, namely: the mosaic model, random sets, models defined by multiple indicators and isofactorial models. In each case, the conditions for the sequential algorithm to honor the model parameters are sought after. Concurrently, the properties of the multivariate distributions are identified and some conceptual impediments are emphasized. In particular, the prior multiple-point statistics are shown to depend on external factors such as the total number of simulated nodes and the number and locations of the samples. As a consequence, common applications such as a flow simulation or a change of support on the realizations may lead to hazardous interpretations.     

Geologic heterogeneity and a comparison of two geostatistical models: Sequential Gaussian and transition probability-based geostatistical simulation

A covariance-based model-fitting approach is often considered valid to represent field spatial variability of hydraulic properties. This study examines the representation of geologic heterogeneity in two types of geostatistical models under the same mean and spatial covariance structure, and subsequently its effect on the hydraulic response to a pumping test based on 3D high-resolution numerical simulation and field data. Two geostatistical simulation methods, sequential Gaussian simulation (SGS) and transition probability indicator simulation (TPROGS) were applied to create conditional realizations of alluvial fan aquifer systems in the Lawrence Livermore National Laboratory (LLNL) area. The simulated K fields were then used in a numerical groundwater flow model to simulate a pumping test performed at the LLNL site. Spatial connectivity measures of high-K materials (channel facies) captured connectivity characteristics of each geostatistical model and revealed that the TPROGS model created an aquifer (channel) network having greater lateral connectivity. SGS realizations neglected important geologic structures associated with channel and overbank (levee) facies, even though the covariance model used to create these realizations provided excellent fits to sample covariances computed from exhaustive samplings of TPROGS realizations. Observed drawdown response in monitoring wells during a pumping test and its numerical simulation shows that in an aquifer system with strongly connected network of high-K materials, the Gaussian approach could not reproduce a similar behavior in simulated drawdown response found in TPROGS case. Overall, the simulated drawdown responses demonstrate significant disagreement between TPROGS and SGS realizations. This study showed that important geologic characteristics may not be captured by a spatial covariance model, even if that model is exhaustively determined and closely fits the exponential function.     

Geologic heterogeneity recognition using discrete wavelet transformation for subsurface flow solute transport simulations

Subsurface flow and solute transport simulations are performed using different scenarios of permeability fields generated from the sequential Gaussian simulation method (SGS), the multiple-point FILTERSIM algorithm and a new multiple-point wavelet-based simulation method (SWS). The SWS method is a multiple-point pattern-based simulation method which uses discrete wavelet transformation for the representation of geologic heterogeneity. For pattern-based simulation, patterns are generated by scanning a training image with a spatial template. The pattern classifications were performed after reducing the dimension of patterns by wavelet decomposition at the suitable scale and by taking only scaling components of wavelet decomposed patterns. The simulation is performed in a sequential manner by finding the best-matched class corresponding to the conditioning data and by randomly sampling a pattern from the best-matched class. The developed method is compared with two other multi-point simulation algorithms, FLTERSIM and SIMPAT. The comparative results revealed that the proposed method is computationally faster than the other two methods while the simulation maps are comparable. Numerical simulations of two flow problems are performed using SGS, SWS and FILTERSIM realizations. The numerical results show a superiority of the SWS method over SGS and FILTERSIM in terms of reproduction of the reference images main features, and agreement with flow and transport results obtained on reference images.     

A multidimensional scaling approach to enforce reproduction of transition probabilities in truncated plurigaussian simulation

Truncated plurigaussian (TPG) simulation is a flexible method for simulating rock types in deposits with complicated ordering structures. The truncation of a multivariate Gaussian distribution controls the proportions and ordering of rock types in the simulation while the variogram for each Gaussian variable controls rock type continuity. The determination of a truncation procedure for complicated geological environments is not trivial. A method for determining the truncation and fitting variograms applicable to any number of rock types and multivariate Gaussian distribution is developed here to address this problem. Multidimensional scaling is applied to place dissimilar categories far apart and similar categories close together. The multivariate space is then mapped using a Voronoi decomposition and rotated to optimize variogram reproduction. A case study simulating geologic layers at a large mineral deposit demonstrates the potential of this method and compares the results with sequential indicator simulation (SIS). Input proportion and transition probability reproduction with TPG is demonstrated to be better than SIS. Variogram reproduction is comparable for both techniques.     

Gamma random field simulation by a covariance matrix transformation method

In studies involving environmental risk assessment, Gaussian random field generators are often used to yield realizations of a Gaussian random field, and then realizations of the non-Gaussian target random field are obtained by an inverse-normal transformation. Such simulation process requires a set of observed data for estimation of the empirical cumulative distribution function (ECDF) and covariance function of the random field under investigation. However, if realizations of a non-Gaussian random field with specific probability density and covariance function are needed, such observed-data-based simulation process will not work when no observed data are available. In this paper we present details of a gamma random field simulation approach which does not require a set of observed data. A key element of the approach lies on the theoretical relationship between the covariance functions of a gamma random field and its corresponding standard normal random field. Through a set of devised simulation scenarios, the proposed technique is shown to be capable of generating realizations of the given gamma random fields.     

Prediction of Remediation of a Heterogeneous Aquifer: A Case Study

Contaminant plumes whose characteristic length is smaller than the horizontal integral scale of the hydraulic conductivity, K, are abundant in shallow, phreatic aquifers. In such cases, the aquifer can be regarded as layered, with K being only a function of the vertical coordinate. The heterogeneity of K has a critical role upon the efficiency of remediation of such sites, for example, by Pump and Treat schemes. The expected efficiency is a random variable, with uncertainty. Quantifying this uncertainty can be of great importance to decision making. In this study, we focus on a case study in the coastal aquifer of Israel and compare two different approaches for constructing realizations of K: continuous and indicator. We observe a significant difference between the constructed realizations, which results in a considerable difference in the predicted remediation efficiency and its uncertainty. Furthermore, we study the effect of conditioning the realizations by a rather limited number of K data points. We find that the conditioning results in a major reduction of the uncertainty. In addition, we compare the results of the transport model to a simplified semi‐analytical solution that is based on assuming radial flow. We find a good agreement with the three‐dimensional numerical model. This result illustrates that the simplified solution can be used for prediction of the remediation efficiency when the flow at the plume vicinity can be regarded as radial.     

Numerical modelling of stream–aquifer interaction: Quantifying the impact of transient streambed permeability and aquifer heterogeneity

Stream–aquifer interaction plays a vital role in the water cycle, and a proper study of this interaction is needed for understanding groundwater recharge, contaminants migration, and for managing surface water and groundwater resources. A model‐based investigation of a field experiment in a riparian zone of the Schwarzbach river, a tributary of the Rhine River in Germany, was conducted to understand stream–aquifer interaction under alternative gaining and losing streamflow conditions. An equivalent streambed permeability, estimated by inverting aquifer responses to flood waves, shows that streambed permeability increased during infiltration of stream water to aquifer and decreased during exfiltration. Aquifer permeability realizations generated by multiple‐point geostatistics exhibit a high degree of heterogeneity and anisotropy. A coupled surface water groundwater flow model was developed incorporating the time‐varying streambed permeability and heterogeneous aquifer permeability realizations. The model was able to reproduce varying pressure heads at two observation wells near the stream over a period of 55 days. A Monte Carlo analysis was also carried out to simulate groundwater flow, its age distribution, and the release of a hypothetical wastewater plume into the aquifer from the stream. Results of this uncertainty analysis suggest (a) stream–aquifer exchange flux during the infiltration periods was constrained by aquifer permeability; (b) during exfiltration, this flux was constrained by the reduced streambed permeability; (c) the effect of temporally variable streambed permeability and aquifer heterogeneity were found important to improve the accurate capture of the uncertainty; and (d) probabilistic infiltration paths in the aquifer reveal that such pathways and the associated prediction of the extent of the contaminant plume are highly dependent on aquifer heterogeneity.     

基于神经网络的随机地震反演方法

针对随机地震反演中存在的两个主要问题,随机实现含有噪声和难以从大量随机实现中挖掘有效信息,提出了一种基于神经网络的随机地震反演方法.通过对多组随机实现及其正演地震数据的计算,构建了基于序贯高斯模拟的训练集.这也为应用神经网络求解地球物理反问题,提供了一种有效建立训练集的方法.较之传统的神经网络反演,这种训练集不仅保证了学习样本具有多样性,同时还引入了空间相关性.数值模拟结果表明,该方法只需要通过单层前馈神经网络,就可以比较有效的解决一个500个阻抗参数的反演问题.     

Testing the correctness of the sequential algorithm for simulating Gaussian random fields

The sequential algorithm is widely used to simulate Gaussian random fields. However, a rigorous application of this algorithm is impractical and some simplifications are required, in particular a moving neighborhood has to be defined. To examine the effect of such restriction on the quality of the realizations, a reference case is presented and several parameters are reviewed, mainly the histogram, variogram, indicator variograms, as well as the ergodic fluctuations in the first and second-order statistics. The study concludes that, even in a favorable case where the simulated domain is large with respect to the range of the model, the realizations may poorly reproduce the second-order statistics and be inconsistent with the stationarity and ergodicity assumptions. Practical tips such as the multiple-grid strategy do not overcome these impediments. Finally, extending the original algorithm by using an ordinary kriging should be avoided, unless an intrinsic random function model is sought after.     

Half-tapering strategy for conditional simulation with large datasets

Gaussian conditional realizations are routinely used for risk assessment and planning in a variety of Earth sciences applications. Assuming a Gaussian random field, conditional realizations can be obtained by first creating unconditional realizations that are then post-conditioned by kriging. Many efficient algorithms are available for the first step, so the bottleneck resides in the second step. Instead of doing the conditional simulations with the desired covariance (F approach) or with a tapered covariance (T approach), we propose to use the taper covariance only in the conditioning step (half-taper or HT approach). This enables to speed up the computations and to reduce memory requirements for the conditioning step but also to keep the right short scale variations in the realizations. A criterion based on mean square error of the simulation is derived to help anticipate the similarity of HT to F. Moreover, an index is used to predict the sparsity of the kriging matrix for the conditioning step. Some guides for the choice of the taper function are discussed. The distributions of a series of 1D, 2D and 3D scalar response functions are compared for F, T and HT approaches. The distributions obtained indicate a much better similarity to F with HT than with T.     

Application of 1D paleo‐fluvial process modelling at a basin scale to augment sparse borehole data: example of a Permian formation in the Galilee Basin,Australia

The heterogeneous hydraulic conductivity (K) in water‐bearing formations controls subsurface flow and solute transport processes. Geostatistical techniques are often employed to characterize the K distribution in space based on the correlation between K measurements. However, at the basin scale, there are often insufficient measurements for inferring the spatial correlation. This is a widespread problem that we address in this study using the example of the Betts Creek Beds (BCB) in the Galilee Basin, Australia. To address the lack of data, we use a 1D stochastic fluvial process‐based model (SFPM) to quantify the total sediment thickness, Z( x ), and the sandstone proportion over the total thickness, P_s( x ), in the BCB. The semivariograms of Z( x ) and P_s( x ) are then extracted and used in sequential Gaussian simulation to construct the 2D spatial distribution of Z( x ) and P_s( x ). P_s( x ) can be converted to a K distribution based on classical averaging methods. The results demonstrate that the combination of SFPM and geostatistical simulation allows for the evaluation of upscaled K distribution with a limited number of K measurements. Copyright © 2015 John Wiley & Sons, Ltd.     

Multiple imputation framework for data assignment in truncated pluri-Gaussian simulation

Truncated pluri-Gaussian simulation (TPGS) is suitable for the simulation of categorical variables that show natural ordering as the TPGS technique can consider transition probabilities. The TPGS assumes that categorical variables are the result of the truncation of underlying latent variables. In practice, only the categorical variables are observed. This translates the practical application of TPGS into a missing data problem in which all latent variables are missing. Latent variables are required at data locations in order to condition categorical realizations to observed categorical data. The imputation of missing latent variables at data locations is often achieved by either assigning constant values or spatially simulating latent variables subject to categorical observations. Realizations of latent variables can be used to condition all model realizations. Using a single realization or a constant value to condition all realizations is the same as assuming that latent variables are known at the data locations and this assumption affects uncertainty near data locations. The techniques for imputation of latent variables in TPGS framework are investigated in this article and their impact on uncertainty of simulated categorical models and possible effects on factors affecting decision making are explored. It is shown that the use of single realization of latent variables leads to underestimation of uncertainty and overestimation of measured resources while the use constant values for latent variables may lead to considerable over or underestimation of measured resources. The results highlight the importance of multiple data imputation in the context of TPGS.     

Building on crossvalidation for increasing the quality of geostatistical modeling

The random function is a mathematical model commonly used in the assessment of uncertainty associated with a spatially correlated attribute that has been partially sampled. There are multiple algorithms for modeling such random functions, all sharing the requirement of specifying various parameters that have critical influence on the results. The importance of finding ways to compare the methods and setting parameters to obtain results that better model uncertainty has increased as these algorithms have grown in number and complexity. Crossvalidation has been used in spatial statistics, mostly in kriging, for the analysis of mean square errors. An appeal of this approach is its ability to work with the same empirical sample available for running the algorithms. This paper goes beyond checking estimates by formulating a function sensitive to conditional bias. Under ideal conditions, such function turns into a straight line, which can be used as a reference for preparing measures of performance. Applied to kriging, deviations from the ideal line provide sensitivity to the semivariogram lacking in crossvalidation of kriging errors and are more sensitive to conditional bias than analyses of errors. In terms of stochastic simulation, in addition to finding better parameters, the deviations allow comparison of the realizations resulting from the applications of different methods. Examples show improvements of about 30% in the deviations and approximately 10% in the square root of mean square errors between reasonable starting modelling and the solutions according to the new criteria.     

Two-phase flow properties of a horizontal fracture: The effect of aperture distribution

A systematic numerical method has been presented to investigate the constitutive relationships between two-phase flow properties of horizontal fractures and aperture distributions. Based on fractal geometry, single rough-walled fractures are generated numerically by modified successive random addition (SRA) method and then aperture distributions with truncated Gaussian distribution are formed by shear displacement between lower and upper surfaces. (The truncated Gaussian distribution is used to describe aperture evolution under different normal stresses.) According to the assumption of two-dimensional porous media and local parallel plate model, invasion percolation approach is employed to model the two-phase flow displacement (imbibition) in generated horizontal fractures, in which capillary forces are dominant over viscous and gravity forces. For truncated Gaussian distributions, constitutive relationships from numerical simulation are compared to closed-form relationships and a good agreement is obtained. The simulation results indicate strong phase interference with the sum of two phase relative permeability values being less than one in the intermediate saturations. It is found that fracture properties related to residual saturations depend on spatial correlation of aperture distributions. Based on the simulation results, we proposed an empirical relationship between the fracture residual-saturation-rated parameters and the corresponding aperture distributions.