Mathematical Geosciences

Geostatistical analyses require an estimation of the covariance structure of a random field and its parameters jointly from noisy data. Whereas in some cases (as in that of a Matérn variogram) a range of structural models can be captured with one or a few parameters, in many other cases it is necessary to consider a discrete set of structural model alternatives, such as drifts and variograms. Ranking these alternatives and identifying the best among them has traditionally been done with the aid of information theoretic or Bayesian model selection criteria. There is an ongoing debate in the literature about the relative merits of these various criteria. We contribute to this discussion by using synthetic data to compare the abilities of two common Bayesian criteria, BIC and KIC, to discriminate between alternative models of drift as a function of sample size when drift and variogram parameters are unknown. Adopting the results of Markov Chain Monte Carlo simulations as reference we confirm that KIC reduces asymptotically to BIC and provides consistently more reliable indications of model quality than does BIC for samples of all sizes. Practical considerations often cause analysts to replace the observed Fisher information matrix entering into KIC with its expected value. Our results show that this causes the performance of KIC to deteriorate with diminishing sample size. These results are equally valid for one and multiple realizations of uncertain data entering into our analysis. Bayesian theory indicates that, in the case of statistically independent and identically distributed data, posterior model probabilities become asymptotically insensitive to prior probabilities as sample size increases. We do not find this to be the case when working with samples taken from an autocorrelated random field.

Geochemical surveys often contain several tens of components, obtained from different horizons and with different analytical techniques. These are used either to obtain elemental concentration maps or to explore links between the variables. The first task involves interpolation, the second task principal component analysis (PCA) or a related technique. Interpolation of all geochemical variables (in wt% or ppm) should guarantee consistent results: At any location, all variables must be positive and sum up to 100 %. This is not ensured by any conventional geostatistical technique. Moreover, the maps should ideally preserve any link present in the data. PCA also presents some problems, derived from the spatial dependence between the observations, and the compositional nature of the data. Log-ratio geostatistical techniques offer a consistent solution to all these problems. Variation-variograms are introduced to capture the spatial dependence structure: These are direct variograms of all possible log ratios of two components. They can be modeled with a function analogous to the linear model of coregionalization (LMC), where for each spatial structure there is an associated variation matrix describing the links between the components. Eigenvalue decompositions of these matrices provide a PCA of that particular spatial scale. The whole data set can then be interpolated by cokriging. Factorial cokriging can also be used to map a certain spatial structure, eventually projected onto those principal components (PCs) of that structure with relevant contribution to the spatial variability. If only one PC is used for a certain structure, the maps obtained represent the spatial variability of a geochemical link between the variables. These procedures and their advantages are illustrated with the horizon C Kola data set, with 25 components and 605 samples covering most of the Kola peninsula (Finland, Norway, Russia).

A novel Gibbs Markov random field for spatial data on Cartesian grids based on the modified planar rotator (MPR) model of statistical physics has been recently introduced for efficient and automatic interpolation of big data sets, such as satellite and radar images. The MPR model does not rely on Gaussian assumptions. Spatial correlations are captured via nearest-neighbor interactions between transformed variables. This allows vectorization of the model which, along with an efficient hybrid Monte Carlo algorithm, leads to fast execution times that scale approximately linearly with system size. The present study takes advantage of the short-range nature of the interactions between the MPR variables to parallelize the algorithm on graphics processing units (GPUs) in the Compute Unified Device Architecture programming environment. It is shown that, for the processors employed, the GPU implementation can lead to impressive computational speedups, up to almost 500 times on large grids, compared to single-processor calculations. Consequently, massive data sets comprising millions of data points can be automatically processed in less than one second on an ordinary GPU.

This study revisits the conditional bias that can be observed with spatial estimators such as kriging. In the geostatistical literature, the term “conditional bias” has been used to describe two different effects: underestimation of high values and overestimation of low values, or the opposite, viz. overestimation of high values and underestimation of low values. To add to the confusion, the smoothing effect of the estimator is always indicated to be the culprit. It seems that geostatisticians have been debating conditional bias since the birth of geostatistics. Is less or more smoothing required to alleviate conditional bias, and which one? This paradox is actually resolved when one considers the different distribution partitions on which conditional expectation can be calculated. Depending on the partitions of the bivariate distribution of true versus estimated values, conditional expectation can be calculated on conditional or marginal distributions. These lead to different types of conditional bias, and smoothing affects them differently. The type based on conditional distributions is smoothing friendly, while the type based on marginal distributions is smoothing adverse. The same estimator can display under- and overestimation, depending on whether a conditional or marginal distribution is considered. It is also observed that all conditional biases, regardless of the bivariate distribution partitions, are greatly affected by the variance of the conditioning data and vary with the sampling. A simple estimator correction can be applied to exactly remove the smoothing-friendly conditional bias in the sample as measured by the slope of the linear regression between the true and estimated values in cross-validation. Over many samplings, it is observed that this cross-validation measure is itself conditionally biased, depending on the variance of the data. On the other hand, the smoothing-adverse type of conditional bias can be corrected by conditional simulation that reproduces the distribution of the data. The results are also biased, depending on the variance of the conditioning data. Correcting for the smoothing-adverse type will worsen the smoothing-friendly type, and vice versa. Both types of conditional bias can be corrected by averaging statistics, or averaging estimates, over multiple samplings.

The mining industry continuously struggles to keep produced tonnages and grades aligned with targets derived from model-based expectations. Deviations often result from the inability to characterise short-term production units accurately based on sparsely distributed exploration data. During operation, the characterisation of short-term production units can be significantly improved when deviations are monitored and integrated back into the underlying grade control model. A previous contribution introduced a novel simulation-based geostatistical approach to repeatedly update the grade control model based on online data from a production monitoring system. The added value of the presented algorithm results from its ability to handle inaccurate observations made on blended material streams originating from two or more extraction points. This contribution further extends previous work studying the relation between system control parameters and algorithm performance. A total of 125 experiments are conducted to quantify the effects of variations in measurement volume, blending ratio and sensor precision. Based on the outcome of the experiments, recommendations are formulated for optimal operation of the monitoring system, guaranteeing the best possible algorithm performance.

Understanding the subsurface structure and function in the near-surface groundwater system, including fluid flow, geomechanical, and weathering processes, requires accurate predictions of the spatial distribution of petrophysical properties, such as rock and fluid (air and water) volumetric fractions. These properties can be predicted from geophysical measurements, such as electrical resistivity tomography and refraction seismic data, by solving a rock physics inverse problem. A Bayesian inversion approach based on a Monte Carlo implementation of the Bayesian update problem is developed to generate multiple realizations of porosity and water saturation conditioned on geophysical data. The model realizations are generated using a geostatistical algorithm and updated according to the ensemble smoother approach, an efficient Bayesian data assimilation technique. The prior distribution includes a spatial correlation function such that the model realizations mimic the geological spatial continuity. The result of the inversion includes a set of realizations of porosity and water saturation, as well as the most likely model and its uncertainty, that are crucial to understand fluid flow, geomechanical, and weathering processes in the critical zone. The proposed approach is validated on two synthetic datasets motivated by the Southern Sierra Critical Zone Observatory and is then applied to data collected on a mountain hillslope near Laramie, Wyoming. The inverted results match the measurements, honor the spatial correlation prior model, and provide geologically realistic petrophysical models of weathered rock at Earth’s surface.

Coal seam degasification and its success are important for controlling methane, and thus for the health and safety of coal miners. During the course of degasification, properties of coal seams change. Thus, the changes in coal reservoir conditions and in-place gas content as well as methane emission potential into mines should be evaluated by examining time-dependent changes and the presence of major heterogeneities and geological discontinuities in the field. In this work, time-lapsed reservoir and fluid storage properties of the New Castle coal seam, Mary Lee/Blue Creek seam, and Jagger seam of Black Warrior Basin, Alabama, were determined from gas and water production history matching and production forecasting of vertical degasification wellbores. These properties were combined with isotherm and other important data to compute gas-in-place (GIP) and its change with time at borehole locations. Time-lapsed training images (TIs) of GIP and GIP difference corresponding to each coal and date were generated by using these point-wise data and Voronoi decomposition on the TI grid, which included faults as discontinuities for expansion of Voronoi regions. Filter-based multiple-point geostatistical simulations, which were preferred in this study due to anisotropies and discontinuities in the area, were used to predict time-lapsed GIP distributions within the study area. Performed simulations were used for mapping spatial time-lapsed methane quantities as well as their uncertainties within the study area. The systematic approach presented in this paper is the first time in literature that history matching, TIs of GIPs and filter simulations are used for degasification performance evaluation and for assessing GIP for mining safety. Results from this study showed that using production history matching of coalbed methane wells to determine time-lapsed reservoir data could be used to compute spatial GIP and representative GIP TIs generated through Voronoi decomposition. Furthermore, performing filter simulations using point-wise data and TIs could be used to predict methane quantity in coal seams subjected to degasification. During the course of the study, it was shown that the material balance of gas produced by wellbores and the GIP reductions in coal seams predicted using filter simulations compared very well, showing the success of filter simulations for continuous variables in this case study. Quantitative results from filter simulations of GIP within the studied area briefly showed that GIP was reduced from an initial ∼73 Bcf (median) to ∼46 Bcf (2011), representing a 37 % decrease and varying spatially through degasification. It is forecasted that there will be an additional ∼2 Bcf reduction in methane quantity between 2011 and 2015. This study and presented results showed that the applied methodology and utilized techniques can be used to map GIP and its change within coal seams after degasification, which can further be used for ventilation design for methane control in coal mines.