Efficient Simulation of (Log)Normal Random Fields for Hydrogeological Applications

Two methods for generating representative realizations from Gaussian and lognormal random field models are studied in this paper, with term representative implying realizations efficiently spanning the range of possible attribute values corresponding to the multivariate (log)normal probability distribution. The first method, already established in the geostatistical literature, is multivariate Latin hypercube sampling, a form of stratified random sampling aiming at marginal stratification of simulated values for each variable involved under the constraint of reproducing a known covariance matrix. The second method, scarcely known in the geostatistical literature, is stratified likelihood sampling, in which representative realizations are generated by exploring in a systematic way the structure of the multivariate distribution function itself. The two sampling methods are employed for generating unconditional realizations of saturated hydraulic conductivity in a hydrogeological context via a synthetic case study involving physically-based simulation of flow and transport in a heterogeneous porous medium; their performance is evaluated for different sample sizes (number of realizations) in terms of the reproduction of ensemble statistics of hydraulic conductivity and solute concentration computed from a very large ensemble set generated via simple random sampling. The results show that both Latin hypercube and stratified likelihood sampling are more efficient than simple random sampling, in that overall they can reproduce to a similar extent statistics of the conductivity and concentration fields, yet with smaller sampling variability than the simple random sampling.     

Modeling spatial uncertainty of heavy metal content in soil by conditional Latin hypercube sampling and geostatistical simulation

This study proposes the method of simulating spatial patterns and quantifying the uncertainty in multivariate distribution of heavy metals (Cr, Cu, Ni, and Zn) by sequential indicator simulation (SIS) combined with conditional Latin hypercube sampling (cLHS) in Changhua County, Taiwan. The cLHS is used for a sampling then for SIS mapping and assessing uncertainties of heavy metal concentrations. The indicator variogram results indicate that the 700 cLHS samples replicate statistical multivariate distribution and spatial structure of the 1,082 samples. Moreover, the SIS realizations based on 700 cLHS samples are more conservative and reliable than those based on 1,082 samples for delineating soil contamination by all heavy metals with the exception of Zn. Given adequate sampling, soil contamination simulation provides sufficient information for delineating contaminated areas and planning environmental management.     

Assessing spatial uncertainty in mapping soil erodibility factor using geostatistical stochastic simulation

Soil erosion is one of most widespread process of degradation. The erodibility of a soil is a measure of its susceptibility to erosion and depends on many soil properties. Soil erodibility factor varies greatly over space and is commonly estimated using the revised universal soil loss equation. Neglecting information about estimation uncertainty may lead to improper decision-making. One geostatistical approach to spatial analysis is sequential Gaussian simulation, which draws alternative, equally probable, joint realizations of a regionalised variable. Differences between the realizations provide a measure of spatial uncertainty and allow us to carry out an error analysis. The objective of this paper was to assess the model output error of soil erodibility resulting from the uncertainties in the input attributes (texture and organic matter). The study area covers about 30 km² (Calabria, southern Italy). Topsoil samples were collected at 175 locations within the study area in 2006 and the main chemical and physical soil properties were determined. As soil textural size fractions are compositional data, the additive-logratio (alr) transformation was used to remove the non-negativity and constant-sum constraints on compositional variables. A Monte Carlo analysis was performed, which consisted of drawing a large number (500) of identically distributed input attributes from the multivariable joint probability distribution function. We incorporated spatial cross-correlation information through joint sequential Gaussian simulation, because model inputs were spatially correlated. The erodibility model was then estimated for each set of the 500 joint realisations of the input variables and the ensemble of the model outputs was used to infer the erodibility probability distribution function. This approach has also allowed for delineating the areas characterised by greater uncertainty and then to suggest efficient supplementary sampling strategies for further improving the precision of K value predictions.     

Integration of Tracer Test Data to Refine Geostatistical Hydraulic Conductivity Fields Using Sequential Self-Calibration Method

On the basis of local measurements of hydraulic conductivity,geostatistical methods have been found to be useful in heterogeneity characterization of a hydraulic conductivity field on a regional scale. However,the methods are not suited to directly integrate dynamic production data,such as,hydraulic head and solute concentration,into the study of conductivity distribution. These data,which record the flow and transport processes in the medium,are closely related to the spatial distribution of hydraulic conductivity. In this study,a three-dimensional gradient-based inverse method-the sequential self-calibration (SSC) method-is developed to calibrate a hydraulic conductivity field,initially generated by a geostatistical simulation method,conditioned on tracer test results. The SSC method can honor both local hydraulic conductivity measurements and tracer test data. The mismatch between the simulated hydraulic conductivity field and the reference true one,measured by its mean square error (MSE),is reduced through the SSC conditional study. In comparison with the unconditional results,the SSC conditional study creates the mean breakthrough curve much closer to the reference true curve,and significantly reduces the prediction uncertainty of the solute transport in the observed locations. Further,the reduction of uncertainty is spatially dependent,which indicates that good locations,geological structure,and boundary conditions will affect the efficiency of the SSC study results.     

An adaptive sampling approach to reduce uncertainty in slope stability analysis

An adaptive sampling approach is proposed, which can sample spatially varying shear strength parameters efficiently to reduce uncertainty in the slope stability analysis. This approach employs a limit equilibrium model and stochastic conditional methodology to determine the likely sampling locations. Karhunen-Loève expansion is used to conduct the conditional Monte Carlo simulation. A first-order analysis is also proposed to ease the computational burden associated with Monte Carlo simulation. These approaches are then tested using borehole data from a field site. Results indicate that the proposed adaptive sampling approach is an effective and efficient sampling scheme for reducing uncertainty in slope stability analysis.     

The use of connectivity clusters in stochastic modelling of highly heterogeneous porous media

Stochastic geostatistical techniques are essential tools for groundwater flow and transport modelling in highly heterogeneous media. Typically, these techniques require massive numbers of realizations to accurately simulate the high variability and account for the uncertainty. These massive numbers of realizations imposed several constraints on the stochastic techniques (e.g. increasing the computational effort, limiting the domain size, grid resolution, time step and convergence issues). Understanding the connectivity of the subsurface layers gives an opportunity to overcome these constraints. This research presents a sampling framework to reduce the number of the required Monte Carlo realizations utilizing the connectivity properties of the hydraulic conductivity distributions in a three-dimensional domain. Different geostatistical distributions were tested in this study including exponential distribution with the Turning Bands (TBM) algorithm and spherical distribution using Sequential Gaussian Simulation (SGSIM). It is found that the total connected fraction of the largest clusters and its tortuosity are highly correlated with the percentage of mass arrival and the first arrival quantiles at different control planes. Applying different sampling techniques together with several indicators suggested that a compact sample representing only 10% of the total number of realizations can be used to produce results that are close to the results of the full set of realizations. Also, the proposed sampling techniques specially utilizing the low conductivity clustering show very promising results in terms of matching the full range of realizations. Finally, the size of selected clusters relative to domain size significantly affects transport characteristics and the connectivity indicators.     

Deterministic and statistical size effect during shearing of granular layer within a micro‐polar hypoplasticity

The paper deals with numerical investigations of a deterministic and statistical size effect in granular bodies during quasi‐static shearing of an infinite layer under plane strain conditions, free dilatancy and constant pressure. For a simulation of the mechanical behaviour of a cohesionless granular material during a monotonous deformation path, a micro‐polar hypoplastic constitutive relation was used which takes into account particle rotations, curvatures, non‐symmetric stresses, couple stresses and the mean grain diameter as a characteristic length. The proposed model captures the essential mechanical features of granular bodies in a wide range of densities and pressures with a single set of constants. In the paper, a deterministic and statistical size effect is analysed. The deterministic calculations were carried out with an uniform distribution of the initial void ratio for four different heights of the granular layer: 5, 50, 500 and 2000 mm. To investigate the statistical size effect, the Monte Carlo method was applied. The random distribution of the initial void ratio was assumed to be spatially correlated. Truncated Gaussian random fields were generated in a granular layer using an original conditional rejection method. The sufficient number of samples was determined by analysing the convergence of the outcomes. In order to reduce the number of realizations without losing the accuracy of the calculations, stratified and Latin hypercube methods were applied. A parametric analysis of these methods was also presented. Some general conclusions were formulated. Copyright © 2007 John Wiley & Sons, Ltd.     

城市暴雨径流模拟的参数不确定性研究

基于暴雨管理模型(SWMM),构建了澳门半岛雅廉访实验小区降雨径流模拟系统,采用generalized likelihood uncertainty estimation(GLUE)方法对复杂降雨径流模型参数的不确定性进行了分析,采用均匀分布为参数的先验分布(prior distribution),选取拉丁超立方法作为采样方法,Nash-Sutcliffe效率系数作为似然评判准则。结果表明该方法可以有效的分析降雨径流模型参数的不确定性,对降雨径流模型参数识别提供了深入分析与理解模型系统的有效途径。     

Wellhead protection area delineation using multiple methods: a case study in Beijing

Determining a wellhead protection area (WHPA) is a basic and important step in protecting groundwater from contamination. Based on availability of data and complexity of hydrogeological conditions, different delineation methods can be adopted. This study’s objective was to examine the relation and difference of the results calculated using six WHPA delineation methods, which ranged from simplified shapes to stochastic models with superior sampling means. All the methods were applied to a well field in Beijing, China. The WHPA probability distribution from a stochastic simulation was used as reference for comparison with other methods. WHPAs calculated using simplified and analytical approaches turned out to be smaller than reference WHPA. Those calculated using the numerical model provided better results but still cannot include a considerable area of reference WHPA. The stochastic approach based on efficient orthogonal Latin hypercube sampling provided an estimation of the probability distribution of WHPA delineation with different degrees of uncertainty, which played an important role in the WHPA delineation of the well field in Beijing. To attain a realistic WHPA delineation in a complex aquifer system, various delineation approaches can be adopted and all results that validate each other must be considered.     

基于LHS方法的地下水随机模拟及开采量风险性 评价

地下水流数值模拟过程中,水文地质参数的不确定性对模拟结果影响很大。以内蒙古鄂尔多斯市某水源地为例,利用拉丁超立方抽样(LHS)方法获得了含水层渗透参数的随机组合,进行地下水流随机模拟。通过对观测资料与计算水位的绝对值平均(MAE)误差和误差均方根(RMSE)的统计分析,获得了模型较为稳定的随机模拟次数是243次。利用该随机模型对水源地设计开采量进行水位预测,并给出允许降深的风险性分布图。结果表明,预测水位和标准差分布符合实际情况,水位降深大于35 m的风险性最大达到75%。     

Representing Spatial Uncertainty Using Distances and Kernels

Assessing uncertainty of a spatial phenomenon requires the analysis of a large number of parameters which must be processed by a transfer function. To capture the possibly of a wide range of uncertainty in the transfer function response, a large set of geostatistical model realizations needs to be processed. Stochastic spatial simulation can rapidly provide multiple, equally probable realizations. However, since the transfer function is often computationally demanding, only a small number of models can be evaluated in practice, and are usually selected through a ranking procedure. Traditional ranking techniques for selection of probabilistic ranges of response (P10, P50 and P90) are highly dependent on the static property used. In this paper, we propose to parameterize the spatial uncertainty represented by a large set of geostatistical realizations through a distance function measuring “dissimilarity” between any two geostatistical realizations. The distance function allows a mapping of the space of uncertainty. The distance can be tailored to the particular problem. The multi-dimensional space of uncertainty can be modeled using kernel techniques, such as kernel principal component analysis (KPCA) or kernel clustering. These tools allow for the selection of a subset of representative realizations containing similar properties to the larger set. Without losing accuracy, decisions and strategies can then be performed applying a transfer function on the subset without the need to exhaustively evaluate each realization. This method is applied to a synthetic oil reservoir, where spatial uncertainty of channel facies is modeled through multiple realizations generated using a multi-point geostatistical algorithm and several training images.     

土壤饱和导水率空间预测的不确定性分析

当土壤转换函数应用于土壤水力性质估计时,对于预测值的不确定性往往容易被忽视。为了有针对性地提出减少这种不确定性的方法和措施,提高土壤转换函数的实际应用能力,以两种现有的土壤转换函数(Vereecken和HYPRES模型)为例,将其应用于山东省平度市土壤饱和导水率的空间预测,并利用拉丁超立方抽样(LHS)方法对预测结果的不确定性进行了分析。结果表明,饱和导水率空间预测的不确定性主要来源于土壤基本性质的空间插值误差和土壤转换函数自身的预测误差。当Vereecken模型应用于饱和导水率空间预测时,预测结果的不确定性主要由土壤基本性质空间插值误差所决定,土壤转换函数预测误差的影响较小,而HYPRES模型则是受二者的双重影响。     

A Fast Approximation for Seismic Inverse Modeling: Adaptive Spatial Resampling

Seismic inverse modeling, which transforms appropriately processed geophysical data into the physical properties of the Earth, is an essential process for reservoir characterization. This paper proposes a work flow based on a Markov chain Monte Carlo method consistent with geology, well-logs, seismic data, and rock-physics information. It uses direct sampling as a multiple-point geostatistical method for generating realizations from the prior distribution, and Metropolis sampling with adaptive spatial resampling to perform an approximate sampling from the posterior distribution, conditioned to the geophysical data. Because it can assess important uncertainties, sampling is a more general approach than just finding the most likely model. However, since rejection sampling requires a large number of evaluations for generating the posterior distribution, it is inefficient and not suitable for reservoir modeling. Metropolis sampling is able to perform an equivalent sampling by forming a Markov chain. The iterative spatial resampling algorithm perturbs realizations of a spatially dependent variable, while preserving its spatial structure by conditioning to subset points. However, in most practical applications, when the subset conditioning points are selected at random, it can get stuck for a very long time in a non-optimal local minimum. In this paper it is demonstrated that adaptive subset sampling improves the efficiency of iterative spatial resampling. Depending on the acceptance/rejection criteria, it is possible to obtain a chain of geostatistical realizations aimed at characterizing the posterior distribution with Metropolis sampling. The validity and applicability of the proposed method are illustrated by results for seismic lithofacies inversion on the Stanford VI synthetic test sets.     

Joint High-Order Simulation of Spatially Correlated Variables Using High-Order Spatial Statistics

Joint geostatistical simulation techniques are used to quantify uncertainty for spatially correlated attributes, including mineral deposits, petroleum reservoirs, hydrogeological horizons, environmental contaminants. Existing joint simulation methods consider only second-order spatial statistics and Gaussian processes. Motivated by the presence of relatively large datasets for multiple correlated variables that typically are available from mineral deposits and the effects of complex spatial connectivity between grades on the subsequent use of simulated realizations, this paper presents a new approach for the joint high-order simulation of spatially correlated random fields. First, a vector random function is orthogonalized with a new decorrelation algorithm into independent factors using the so-termed diagonal domination condition of high-order cumulants. Each of the factors is then simulated independently using a high-order univariate simulation method on the basis of high-order spatial cumulants and Legendre polynomials. Finally, attributes of interest are reconstructed through the back-transformation of the simulated factors. In contrast to state-of-the-art methods, the decorrelation step of the proposed approach not only considers the covariance matrix, but also high-order statistics to obtain independent non-Gaussian factors. The intricacies of the application of the proposed method are shown with a dataset from a multi-element iron ore deposit. The application shows the reproduction of high-order spatial statistics of available data by the jointly simulated attributes.     

A flow-based pattern recognition algorithm for rapid quantification of geologic uncertainty

Geologic uncertainties and limited well data often render recovery forecasting a difficult undertaking in typical appraisal and early development settings. Recent advances in geologic modeling algorithms permit automation of the model generation process via macros and geostatistical tools. This allows rapid construction of multiple alternative geologic realizations. Despite the advances in geologic modeling, computation of the reservoir dynamic response via full-physics reservoir simulation remains a computationally expensive task. Therefore, only a few of the many probable realizations are simulated in practice. Experimental design techniques typically focus on a few discrete geologic realizations as they are inherently more suitable for continuous engineering parameters and can only crudely approximate the impact of geology. A flow-based pattern recognition algorithm (FPRA) has been developed for quantifying the forecast uncertainty as an alternative. The proposed algorithm relies on the rapid characterization of the geologic uncertainty space represented by an ensemble of sufficiently diverse static model realizations. FPRA characterizes the geologic uncertainty space by calculating connectivity distances, which quantify how different each individual realization is from all others in terms of recovery response. Fast streamline simulations are employed in evaluating these distances. By applying pattern recognition techniques to connectivity distances, a few representative realizations are identified within the model ensemble for full-physics simulation. In turn, the recovery factor probability distribution is derived from these intelligently selected simulation runs. Here, FPRA is tested on an example case where the objective is to accurately compute the recovery factor statistics as a function of geologic uncertainty in a channelized turbidite reservoir. Recovery factor cumulative distribution functions computed by FPRA compare well to the one computed via exhaustive full-physics simulations.     

基于地质统计学的NDVI图像估值技术

将疏采样后的NDVI图像作为未受云层影响的已知数据,分别用普通克里格、泛克里格、指示克里格和序贯指示模拟对NDVI图像进行恢复并比较其效果。研究发现,各种克里格法对NDVI图像的估值效果由高到低依次为泛克里格、普通克里格、指示克里格,通常计算方便的普通克里格法就能够满足图像恢复所要求的精度;普通克里格方差和泛克里格方差只能反映数据的构型,不能很好地衡量估值图像的不确定性,指示克里格的条件方差的分布和实际误差的分布基本一致,能够较好地衡量估值图像的不确定性,并且其大小与NDVI影像数据的不确定性大小的分布一致。序贯指示模拟得到的多个等概率实现表现出很大的空间变异性,多个实现的均值图像光滑效应明显,估值精度不高,但是多个实现的方差分布可以很好地表征空间数据的不确定性分布。     

Application of multiple-point geostatistics on modelling groundwater flow and transport in a cross-bedded aquifer (Belgium)

Sedimentological processes often result in complex three-dimensional subsurface heterogeneity of hydrogeological parameter values. Variogram-based stochastic approaches are often not able to describe heterogeneity in such complex geological environments. This work shows how multiple-point geostatistics can be applied in a realistic hydrogeological application to determine the impact of complex geological heterogeneity on groundwater flow and transport. The approach is applied to a real aquifer in Belgium that exhibits a complex sedimentary heterogeneity and anisotropy. A training image is constructed based on geological and hydrogeological field data. Multiple-point statistics are borrowed from this training image to simulate hydrofacies occurrence, while intrafacies permeability variability is simulated using conventional variogram-based geostatistical methods. The simulated hydraulic conductivity realizations are used as input to a groundwater flow and transport model to investigate the effect of small-scale sedimentary heterogeneity on contaminant plume migration. Results show that small-scale sedimentary heterogeneity has a significant effect on contaminant transport in the studied aquifer. The uncertainty on the spatial facies distribution and intrafacies hydraulic conductivity distribution results in a significant uncertainty on the calculated concentration distribution. Comparison with standard variogram-based techniques shows that multiple-point geostatistics allow better reproduction of irregularly shaped low-permeability clay drapes that influence solute transport.     

Comparison of three geostatistical methods for hydrofacies simulation: a test on alluvial sediments

The hydrodispersive properties of porous sediments are strongly influenced by the heterogeneity at fine scales, which can be modeled by geostatistical simulations. In order to improve the assessment of the properties of three different geostatistical simulation methods (Sequential indicator simulation, SISIM; Transition probability geostatistical simulation, T-PROGS; Multiple point simulation, MPS) a comparison test at different scales was performed for a well-exposed aquifer analogue. In the analysed volume (approximately 30,000?m³) four operative hydrofacies have been recognised: very fine sand and silt, sand, gravelly sand and open framework gravel. Several equiprobable realizations were computed with SISIM, MPS and T-PROGS for a test volume of approximately 400?m³ and for the entire volume, and the different outcomes were compared with visual inspection and connectivity analysis of the very or poorly permeable structures. The comparison of the different simulations shows that the geological model is best reproduced when the simulations are realised separately for each highest rank depositional element and subsequently merged. Moreover, the three methods yield different images of the volume; in particular MPS is efficient in mapping the geometries of the most represented hydrofacies, whereas SISIM and T-PROGS can account for the distribution of the less represented facies.     

DNAPLs污染含水层多相流数值模拟模型的替代模型

针对表面活性剂强化的重非水相流体(DNAPLs)污染的含水层修复问题,在建立多相流数值模拟模型的基础上,应用拉丁超立方采样(LHS)方法,在多相流模拟模型可控输入变量的可行域内采样,有效提高了采样效率和覆盖程度。根据采集的样品数据集,运用多元回归分析方法建立多相流模拟模型的替代模型--双响应面模型,为DNAPLs污染含水层修复过程的优化设计的耦合技术探索新的理论和方法。经检验,替代模型计算结果的相对误差均小于10%,精度较高,说明其在功能上充分逼近模拟模型。运用替代模型实现模拟模型与优化模型的连接,可以大幅度减少优化模型计算过程中直接多次反复调用模拟模型所引起的庞大计算负荷。     

Robust quantification of parametric uncertainty for surfactant–polymer flooding

Uncertainty in surfactant–polymer flooding is an important challenge to the wide-scale implementation of this process. Any successful design of this enhanced oil recovery process will necessitate a good understanding of uncertainty. Thus, it is essential to have the ability to quantify this uncertainty in an efficient manner. Monte Carlo simulation is the traditional uncertainty quantification approach that is used for quantifying parametric uncertainty. However, the convergence of Monte Carlo simulation is relatively low, requiring a large number of realizations to converge. This study proposes the use of the probabilistic collocation method in parametric uncertainty quantification for surfactant–polymer flooding using four synthetic reservoir models. Four sources of uncertainty were considered: the chemical flood residual oil saturation, surfactant and polymer adsorption, and the polymer viscosity multiplier. The output parameter approximated is the recovery factor. The output metrics were the input–output model response relationship, the probability density function, and the first two moments. These were compared with the results obtained from Monte Carlo simulation over a large number of realizations. Two methods for solving for the coefficients of the output parameter polynomial chaos expansion are compared: Gaussian quadrature and linear regression. The linear regression approach used two types of sampling: full-tensor product nodes and Chebyshev-derived nodes. In general, the probabilistic collocation method was applied successfully to quantify the uncertainty in the recovery factor. Applying the method using the Gaussian quadrature produced more accurate results compared with using the linear regression with full-tensor product nodes. Applying the method using the linear regression with Chebyshev derived sampling also performed relatively well. Possible enhancements to improve the performance of the probabilistic collocation method were discussed. These enhancements include improved sparse sampling, approximation order-independent sampling, and using arbitrary random input distribution that could be more representative of reality.