Search

Mineral resources are typically quantified by estimating the grade–tonnage curve for different resource categories. Statutory resource assessment reports (e.g., NI-43-101), however, do not include a measure of uncertainty for the disclosed resources despite the major investments required for mining projects. Although conditional simulation can provide confidence intervals (CIs) for resource estimation, it requires a strong stationarity assumption and depends on the variogram model selected, which is often poorly defined with the available data. In order to avoid these limitations, this research proposes the use and comparison of two machine learning (ML) methods, multiple linear regression and a multilayer neural network, to generate tonnage curves and their CIs directly from the data. The classical variogram modeling step is replaced by the specification of intervals for each parameter of the selected variogram model. The learning is carried out in a perfectly controlled environment using simulations with known variograms. Numerous reference deposits are sampled, and for each one, a series of conditional realizations define the mean tonnage and CI curves. Different statistics computed for the entire data set are used as input to predict the tonnage and CI curves by the ML methods. The results indicate that there are no significant differences between the ML methods. In addition, ML resource predictions outperform those obtained with ordinary kriging, constrained kriging, uniform conditioning and indirect lognormal correction, being surpassed only by the discrete Gaussian model. Nevertheless, these predictors were favored by the use of true variogram models. Moreover, the coverage probabilities of different CIs reach the nominal levels indicating adequate resource uncertainty quantification. Finally, two case studies validate the effectiveness of the proposed approach for tonnage prediction and uncertainty quantification.

The superposition principle is widely used to accelerate the computation of hourly fluid temperatures in ground heat exchangers. This work presents a first attempt at adapting this technique to non-stationary situations of ground-source heat pump operations that involve time-variant parameters, such as flow rates. An advanced response model is first employed to evaluate the outlet temperature signals corresponding to different sets of constant operating conditions. The transitions are then considered by combining the functions linearly and by using a scaled correction function. Eight synthetic validation studies were conducted, representing the dynamic operation of a closed-loop and a standing column well with time-variant circulation, pumping, and discharge rates, and groundwater velocities. It is shown that the developed method reproduces reference numerical results with mean absolute and relative errors that are lower than 0.04 ° C and 0.68%, respectively, and achieves a 30-day simulation with a 5-min time step in a few hundredths of a second. This represents a reduction of five to seven orders of magnitude in computing times and demonstrates the potential of the proposed method to accelerate the simulation of ground heat exchangers that operate with time-variant flow rates.

Obtaining representative samples is of the uttermost importance for modelling mineral deposits. On one hand, Gy’s formula makes it possible to compute relative error variance for a specified target grade on the stringent assumption of perfect homogenization. On the other hand, duplicates (field, rejects and pulp) are now routinely used in QA/QC. No direct links have been established between Gy’s relative variance and the main duplicates statistics used in QA/QC. This paper presents a simple multiplicative error model that makes it possible to predict from Gy’s formula the standard deviation of duplicates differences and duplicates correlation for any distribution. Assuming lognormality, a test is derived to predict whether the sample preparation procedure is likely to meet the target proportion of duplicates with half absolute relative difference (HARD) of less than 10%. The test and the predicted statistics can be used to predict duplicates’ QA/QC results and make necessary adjustments early, before any QA/QC data is acquired. The relationships are tested and validated by simulations and a series of deposit data extracted from NI43-101 published reports.

The Sequential Spectral Turning Bands Method (S-STBM) builds Gaussian random fields (GRF) calibrated to desired response functions. An interesting application of S-STBM concerns the simulation of GRF subject to inequality constraints. S-STBM works by choosing the phase of each cosine function of the STBM algorithm instead of perturbating nodes of the GRF many thousand times using conditional distributions as in Gibbs sampler. Each chosen phase increasingly constrains the nodes to the desired inequalities. A method based on the sequential Gaussian simulation is introduced to accelerate convergence at the end of the process. Examples shown compare S-STBM approach to Gibbs sampler. Orders of magnitude reduction in computation time is achieved with our spectral method. Furthermore, examples show that the phase selection has no significant influence on the spatial correlation. Our approach is easily generalized to pluriGaussian simulations. Compared to Gibbs sampler, S-STBM is not limited to small systems (no memory limitation) and its complexity of O(n) makes it an efficient tool to simulate large GRF subject to inequality constraints.

In advanced exploration projects or operating mines, the process of allocating capital for infill drilling programs is a significant and recurrent challenge. Within a large company, the different mine sites and projects compete for the available funds for drilling. To maximize a project’s value to its company, a drillhole location optimizer can be used as an objective tool to compare drilling campaigns. The fast semi-greedy optimizer presented here can allow for the obtention of close to optimal solutions to the coverage problem with up to three orders of magnitude less computing time needed than with integer programming. The heuristic approach is flexible as it allows dynamic updating of block values once new drillholes are selected in the solution, as opposed to existing methods based on static block values. The block values used for optimization incorporate kriging estimate and variance, estimate of indicator at cutoff grade and distances to existing or newly selected drillholes. The heuristic approach tends to locate new drillholes within the maximum risk areas, i.e., within less informed zones predicted as being ore zones. Applied to different deposits, it enables, after suitable normalization, comparison of different drilling campaigns and allocation of budgets accordingly.

In most NI-43-101 resource assessment reports the prediction of global in situ resources is performed by either inverse distance weighting, ordinary kriging (OK) or uniform conditioning (UC). These methods have known drawbacks: OK estimates are oversmoothed, and UC necessitates an additional step to localize resources within panels. An alternative, named constrained kriging (CK), enables to circumvent the smoothing issue of OK by imposing the desired theoretical variance to the interpolated variable. CK is not used in NI-43-101 reports, possibly due to a lack of real application examples and little detailed study of its properties. This paper seeks to fill the gap by comparing the prediction performance for global resources of OK, UC and CK on a synthetic lognormal dataset and two real datasets, the Walker Lake and a gold deposit. Results indicate that CK, although being slightly less precise than OK, provides better predictions of grade-tonnage curves than OK and predictions comparable to UC, a remarkable achievement considering that UC is a widespread nonlinear method specifically designed to predict recovery functions. CK is also shown to provide resource estimates more robust than UC with respect to the variogram model specification. Hence, CK appears as a valuable tool allowing simultaneously to localize resources and easily account for change of support in resources estimation.

Significant energy savings can be generated by standing column wells. Despite a growing interest in these ground heat exchangers, very few simulation models have been validated against reliable fie...

The analysis of thermal response tests by optimization-based methods can lead to highly correlated parameters and to a relatively uncertain identification of the thermal parameters. The multi-objective optimization strategy examined here combines the common root mean square error to describe the misfit between the measured and simulated fluid temperature with the root mean square error between the temperature derivative of the same variables. The final objective function combines also the partial correlation coefficients of the temperature and derivative signals. The four objective functions are combined by scalarizing the multi-objective problem. The method proposed here was tested using three real data sets produced by the thermal response tests performed on a helical shaped pipe and on two single U-bent ground heat exchangers. The multi-objective strategy did not lead to different thermal parameters for the three tests, but it did generate a significant reduction in the correlation and uncertainty of the parameters. In addition, a more compact optimum shape can lead to a more rapid convergence time.

This short note discusses the admissibility of a bivariate correlation model recently proposed in this journal. It shows both theoretically and with counter-examples that the proposed model is not admissible in a variety of situations, so that it should be used with caution. Closed form expression of admissibility condition for the bivariate exponential correlation model is derived.

Ground source heat pump systems are among the most energy-efficient heating and cooling technologies. Their performance is strongly related to the accuracy of the ground heat exchanger sizing, hence requiring the forecast of the system’s temperature evolution in response to the anticipated thermal loads. Through this process, simulation techniques that make use of the superposition principle are commonly used to reduce the computational burden. In their current state, these techniques are however only suitable for addressing linear and stationary problems and do not apply to fundamental non stationary situations related to ground source heat pumps operation that involve time-variant parameters.The present work addresses this issue by introducing a novel method based on the principle of superposition that tackles the fast evaluation of the temperature of a closed-loop ground heat exchanger operating with a dynamic heat load as well as time-variant circulation flow rates. The developed method relies on the non stationary combination, a technique borrowed from the field of seismic data processing. This technique achieves discontinuous transitions of convolution products that can be smoothened near transition times by realizing a linear interpolation over the duration of the fluid residence time.The accuracy and efficiency of the proposed method are verified by comparing its results with those provided by reference 3D finite-elements models developed in the Comsol Multiphysics environment. For this purpose, comparative simulations representing the non stationary operation of a closed-loop system having time-variant circulation flow rates are conducted. The case of a single well is first investigated, followed by a borefield of eight wells to demonstrate the validity of the method in both scenarios.Findings indicate that the proposed method can reproduce the reference results with a mean absolute error that is lower than 0.02 °C, and that it is faster than the numerical models by several orders of magnitude. These findings suggest that a broader range of operating scenarios can be handled by highly efficient simulation tools based on the superposition principle, which could foster the development of optimal operating strategies and lead to enhanced overall performances of ground source heat pump systems.

Standing column wells (SCWs) are efficient ground heat exchangers that use local groundwater as a heat source/sink for heating and cooling buildings. In a SCW, high heat exchange rates are achieved by recirculating groundwater in a single deep (75 m to 450 m) and uncased borehole. Discharging (“bleeding”) a small amount of the pumped water outside the SCW also allows maintaining the groundwater temperature within the heat pump’s operational range during peak demand periods. This strategy has been identified as the most significant parameter of SCW operation and is associated with reductions in total length, surface and cost of the borehole heat exchanger compared with the more common closed-loop systems.This work aims at improving knowledge of the dynamic mass and heat transfer processes involved in SCW operation, in order to promote adoption of this energy-efficient technology and encourage good practice. To this end, data is collected using an experimental SCW system located near the city of Montreal, Canada, and made of a 215-m-deep SCW and a 150-m-deep injection well available for discharge of bleed water. The wells are also connected to a large-scale geothermal laboratory designed and equipped to mimic the heating and cooling operation of a small commercial building. First, an advanced finite-element model coupling advection-diffusion of heat and groundwater flow within a SCW and the surrounding ground is developed in the Comsol Multiphysics environment and is validated using experimental datasets collected through downhole temperature measurements, a pumping test, a thermal response test as well as 25 days of winter operation. The numerical model is then used to evaluate the impact of the pumping arrangement and bedrock fracturation on the well’s outlet temperature. Secondly, the operational parameters logged during the dynamic heat extraction test are analysed to provide insight about various operating strategies and their effect on the system’s performance.The work conducted so far demonstrates that the proposed finite-element model reproduces the hydraulic and thermal behaviours of a SCW with satisfying accuracy. Numerical results suggest that placing the submersible pump near the top of the well avoids installation and maintenance difficulties without compromising heat pump operation compared with the usual reverse configuration. It is also shown that deep fractured zones are beneficial to heat pump operation in heating mode, whereas near-surface fracturing tends to impair the performance of the system throughout winter as it eventually favours recharge of the well with colder water. At last, analysis of the winter test data indicates the effectiveness of a three-level bleed control and on-off sequence for maintaining the groundwater temperature above the freezing point, while minimizing the volume of discharged water and allowing to reach a 160 W/m heat extraction rate.

This paper proposes a new approach to the mining exploration drillholes positioning problem (DPP) that incorporates both geostatistical and optimization techniques. A metaheuristic was developed to solve the DPP taking into account an uncertainty index that quantifies the reliability of the current interpretation of the mineral deposit. The uncertainty index was calculated from multiple deposit realizations obtained by truncated Gaussian simulations conditional to the available drillholes samplings. A linear programming model was defined to select the subset of future drillholes that maximizes coverage of the uncertainty. A Tabu Search algorithm was developed to solve large instances of this set partitioning problem. The proposed Tabu Search algorithm is shown to provide good quality solutions approaching 95% of the optimal solution in a reasonable computing time, allowing close to optimal coverage of uncertainty for a fixed investment in drilling.

A dynamic thermal response test is conducted on an experimental standing column well located in the Nicolet Formation to evaluate the impact of bleed rate of 7.5 L/min increases the apparent thermal conductivity by 2.5 times compared with a no-bleed scenario.

A compact block matrix formulation allowing fast construction of the g-function of a ground heat exchanger is presented. This new formulation is nor sequential nor iterative, doesn't require use of Laplace or Fourier transform and provides a g-function by solving a single system of linear equations assembled in a block matrix system. A method to accelerate the integration of the finite-line source model based on Chebyshev polynomials is also introduced. Although it suffers from a fixed cost in terms of computational time, this approach allows to speed up the g-function assessment even further when used jointly with the block matrix formulation on large fields. By using both strategies, constructing the g-function of a ground heat exchanger composed of 50 regularly spaced boreholes can be achieved in less than half-second while committing only a small relative error. The speed and compactness of the block matrix formulation could be useful to design ground heat exchangers with optimization-based algorithms, which can require the assessment of several thousand g-functions.

Methods to simulate facies (or categorical) fields are numerous. However, calibration of simulated facies fields to large-scale or dynamic data still remains an important challenge due to the discrete nature of the fields, the non-linearity of the response with respect to the facies fields, and the non-derivability of the objective function used in calibration. A new gradual deformation method (GDM) is presented and tested for the calibration of facies realizations obtained by patch-multipoint simulation (MPS). The proposed method borrows ideas from pluriGaussian simulation, evolutionary algorithms, and GDM. Various test cases are considered: proportion maps, section of seismic amplitudes, inlet to outlet travel time along the shortest path, and water-cut curves obtained with a flow simulator. Both conditional/unconditional MPS simulations and 2D/3D problems are considered. In all studied test cases, the new GDM approach has provided excellent calibration to the target variables. The method is general as it can be used in conjunction with any facies simulator.

Understanding the geological uncertainty of hydrostratigraphic models is important for risk assessment in hydrogeology. An important feature of sedimentary deposits is the directional ordering of hydrostratigraphic units (HSU). Geostatistical simulation methods propose efficient algorithm for assessing HSU uncertainty. Among different geostatistical methods to simulate categorical data, Bayesian maximum entropy method (BME) and its simplified version Markov-type categorical prediction (MCP) present interesting features. In particular, the zero-forcing property of BME and MCP can provide a valuable constrain on directional properties. We illustrate the ability of MCP to simulate vertically ordered units. A regional hydrostratigraphic system with 11 HSU and different abundances is used. The transitional deterministic model of this system presents lateral variations and vertical ordering. The set of 66 (11 × 12/2) bivariate probability functions is directly calculated on the deterministic model with fast Fourier transform. Despite the trends present in the deterministic model, MCP is unbiased for the HSU proportions in the non-conditional case. In the conditional cases, MCP proved robust to datasets over-representing some HSU. The inter-realizations variability is shown to closely follow the amount and quality of data provided. Our results with different conditioning datasets show that MCP replicates adequately the directional units arrangement. Thus, MCP appears to be a practical method for generating stochastic models in a 3D hydrostratigraphic context.

Fracture spacing (s) distribution and its dependence on bed thickness (t) are investigated for opening-mode joints in sandstone interlayered with shale in the Sainte-Anne des Monts area (the Gaspe Peninsula), Quebec, Canada. The s values of cross joints, which extend across intervals between long parallel systematic joints and abut these joints at about 90° angles, are consistently larger than those of systematic joints by a mean factor of 1.81. The s-t relationship follows a power law: s = mt1-1/k, where m equals half the ratio of layer tensile strength to interface shear stress, and k is the Weibull modulus. The m value is significantly larger for cross joints than systematic joints, indicating that the interfacial shear stress was direction-dependent during the deformation. The Weibull modulus of sandstone was larger during the formation of systematic joints than that of cross joints by a factor of ~1.2, indicating that the number and size of flaws in a unit of volume of the sandstone increased during progressive deformation. Furthermore, the compilation of s-t data from the literature yields: m = 2.53 ± 2.12 and k = 3.09 ± 1.24 for sandstones, and m = 3.53 ± 3.11 and k = 3.54 ± 1.33 for carbonates. The analysis also raises doubt about the existence of real joint saturation in sedimentary rocks.

A 3D finite element model of circular borehole cluster connected in series is used as a response model to generate normalized transfer functions of a borehole thermal energy storage system under the influence of ground surface temperature variations and groundwater flow. Two response functions are obtained by convolving transfer functions of the system's outlet fluid temperature with two incremental input functions describing the temperature variation of the inlet and outlet fluid and variations of the ground surface temperature. Superposition is applied to obtain the resulting fluid temperature with respect to thermal inputs at the inlet fluid and the ground surface over the course of ten years and results are compared for various groundwater velocities. This work demonstrates that the combined effect of groundwater flow and ambient air temperature variations can significantly decrease the performance of a BTES system. The methodology used can be extended to simulate complex system.

In mineral resource evaluations, geostatistical methods known as global change of support allow prediction of the theoretical histogram and gradetonnage curves prior to interpolations or simulations of grades. Two methods commonly used by professionals to guide the choice of interpolation parameters and assess results are the discrete Gaussian model (DGM) and the indirect lognormal correction (IndLog). These models rely upon an estimate of the dispersion variance of the blocks, which is derived by numerical integration of the variogram model over a discretized block. Due to difficulties in obtaining well-formed traditional experimental variograms (especially in the presence of outliers and limited clustered data), many professionals prefer to use 'normalized' variograms such as correlogram (non-ergodic variogram), pairwise-relative, or variogram of the normal score transform. A series of simulations with different grade distributions and variogram models are used to assess the performances and robustness of the various variogram estimators with respect to the DGM and IndLog global change of support. Our results show that the traditional variogram, the correlogram, and normal score variogram have better performances, compared to pairwise, for both DGM and IndLog. Moreover, DGM provided better results than IndLog for the grade distributions that are not strictly lognormal. These findings provide valuable guides for geostatistics practitioners.

Quasi-point hydraulic properties (K) measured locally under laboratory or field conditions need to be upscaled to block-scale K-tensors for use in flow simulators. The upscaled model also needs to be calibrated to hydraulic head observations. The calibration must preserve spatial covariance, cross-covariance and non-linear relations between tensor components. We apply a new upscaling method that allows to compute and model the covariance between block K-tensor components. We use a gradual deformation method for calibration of simulated K-tensor fields to measured head data. Our method incorporates a new bivariate transform that preserves the non-linear relations between K-tensor components. The ensemble of calibrated realizations allows quantification of uncertainty of groundwater flow models. A comparison with PEST on a test case defining capture zones for water supply wells shows that our method calibrates better to measured heads than PEST, provides more realistic K-tensors and results in larger capture zones.

A new approach is described to allow conditioning to both hard data (HD) and soft data for a patch- and distance-based multiple-point geostatistical simulation. The multinomial logistic regression is used to quantify the link between HD and soft data. The soft data is converted by the logistic regression classifier into as many probability fields as there are categories. The local category proportions are used and compared to the average category probabilities within the patch. The conditioning to HD is obtained using alternative training images and by imposing large relative weights to HD. The conditioning to soft data is obtained by measuring the probability–proportion patch distance. Both 2D and 3D cases are considered. Synthetic cases show that a stationary TI can generate non-stationary realizations reproducing the HD, keeping the texture indicated by the TI and following the trends identified in probability maps obtained from soft data. A real case study, the Mallik methane-hydrate field, shows perfect reproduction of HD while keeping a good reproduction of the TI texture and probability trends.

A new algorithm for the optimal stope design problem is proposed. It is based on a previous methodology developed by Bai et al. (Comput Geosci 52:361–371, 2013a) where a cylindrical coordinate system is used to define geomechanical restrictions and to find the optimal stope around an initial raise. The new algorithm extends this work using an integer programming formulation and a new set of constraints, aimed to solve geomechanical issues present on the original methodology. The new formulation is tested on two synthetic and one real deposits. An economic, geomechanical and feasibility analysis is performed, comparing the new results with Bai’s methodology. This methodology achieves better stope designs in terms of geomechanical stability and wall regularity, generating feasible stopes for real use. It also allows further extensions to incorporate other geometrical constraints in order to obtain more regular stope designs.

A fast and efficient method to simulate multivariate fields with non-linear models of coregionalization (N-LMC) is described. The method generalizes FFTMA to the multivariate simulation of the N-LMC with symmetric cross-covariances, hence the name GFFTMA. It allows us for example to use an exponential model as the direct covariance for the main variable, a Cauchy model for the secondary variable and a K-Bessel model for the cross-covariance. Each covariance and cross-covariance are Fast Fourier Transformed (FFT) to get the discrete spectral densities. Then the spectral matrix is eigen-decomposed at each frequency separately to provide the square root matrix and to enforce positive-definiteness in cases where small negative eigenvalues are found. Finally the simulated spectrum is obtained as multiplication of the root matrix and the white noise coefficients. The method is particularly fast for covariances having derivatives at the origin and/or for covariances with long range. Hence, two-variables' 2D fields of 100 million pixels with all-Gaussian or all-cubic covariances and cross-covariance are both simulated in less than 200s. The CPU-time increases only as Nlog(N) (N, the number of points to simulate). Additional realizations are obtained at a low marginal cost as the eigen-decomposition step needs to be done only once for the first realization. The main limitation of the approach is its rather stringent memory requirement. Synthetic examples illustrate the simulations of N-LMC with two and three variables for different combinations of the seven available models. It shows that the theoretical models are all well reproduced. An illustrative case-study on overburden thickness simulation is provided where the secondary information consists of a latent Gaussian variable identifying the geological domain. HighlightsA FFT based method is proposed to simulate non-linear coregionalizations (N-LMC).The method uses FFT and eigen-decomposition of spectral matrices at each frequency.Realizations fit very well the target N-LCM.Computing time is very fast and shows a Nlog(N) relation with size N.An approach based on Cholesky is suggested to get simulated values at sample points.

The remediation of contaminated sites is recognized as a means towards the sustainable development of an expanding urban space on the global scale. However, it is often hindered by financial and environmental risks caused by undetected contamination or unreliable site characterization. Uncertainty in site characterization studies is suspected to be a cause of this problem. Therefore, there is a need for formally defining and quantifying the uncertainty linked to estimates of the spatial distribution of contamination and the volume of contaminated soils at polluted sites. The objective of this study was to develop a geostatistical procedure for the 3-D analysis of a multi-element contamination in order to: 1) optimize the location and the number of sampling stations during a complementary sampling campaign, 2) estimate the volume of contaminated soil for a given data set and 3) quantify the uncertainty of this estimate. This procedure is innovative because it takes into account, by a transformation in pollution indices, the presence of a multi-element contamination and applies to a set of 3-D data. Furthermore, this procedure offers a formal framework that adapts the decisional bounds to the study site to help the decision-maker. The proposed procedure was applied to an urban site characterized by the presence of a strongly heterogeneous fill. From a subset initially constituted of 8 sampling stations, supplemental sampling stations were selected within a complementary data set and added one by one to the iterative geostatistical simulations. The completeness of the sampling plan was determined upon reaching stability in the relative differences of both global uncertainty and global variance between two successive steps of the optimization. According to the procedure presented, a sampling plan could be considered optimized and complete with 19 sampling stations for this study site (comparatively to 38 sampling stations for the total set obtained by systematic sampling). For 19 sampling stations, the global uncertainty was determined to be 20% (comparatively to 16% for the total set).

A coupled model of a standing column well is developed to evaluate the influence of groundwater flow in fractured aquifers. Heat transfer and groundwater flow within a standing column well and its surrounding ground is modeled by means of a resistance and capacity network. The work demonstrates that the presence of a single fracture zone embedded within the bedrock has a significant beneficial effect on the performance of standing column wells. Temperature differences between the homogeneous and fractured aquifer were particularly important when the well was operated at typical bleed ratios. The results show that the model developed in this paper is in good agreement with numerical reference solutions.

Bayesian inference has tremendous potential for thermal response test analysis, as it provides uncertainty metrics that are useful for the design of ground-coupled heat pump systems. The inference process is computationally heavy and has so far been limited to a few thermal parameters and under the unrealistic assumption of residuals’ independence. In this work, a new closed-form expression of the likelihood and an improved artificial neural network are used to speed up Bayesian inference and consider the strong temporal correlation of the residuals. This efficient strategy allowed the robust inference of the joint distribution of five parameters. Using data measured during a real test of 168 h, this work shows that it is possible to robustly identify the volumetric heat capacity of the ground and grout with an uncertainty of 16.3 and 13.8%, a significant improvement. For the specific data used, it is shown that with independence assumption, some parameters are clearly unrealistic, a problem not encountered when the correlation of the residuals is considered. The impact of the interpretation model, of the test duration and of the sampling frequency was also assessed and illustrated by the sizing of a ground heat exchanger. Results reveal that joint identification of some thermal parameters cannot be achieved reliably by the finite line source model, that duration of thermal response tests should be at least 72 h to avoid large uncertainties on the parameters, and that recording temperature every 2 min degrades the identification of the volumetric heat capacity.

The design of pushbacks in an open pit mine has a significant impact on the mine's profitability. Automatic generation of practical pushbacks is a highly desirable feature, but current automatic solutions fail to sufficiently account for complex geometric requirements of pushbacks, including slopes, phase bench and bottom width, smoothness, and continuity. In this paper, we present a tool to fill this gap. Our proposed algorithm is based on modification of sets of blocks obtained by parametric optimization of the pit using a maximum flow method. A set of geometric operators is developed to modify the sets to present a feasible geometry for mining. The geometric operators are essentially derived from mathematical morphological tools. Case studies show that the proposed method generates practical pushback designs that meet all geometric constraints. The algorithm can be used to create a solution for medium-size pits in minutes. This significantly improves the efficiency of designing pushbacks for open pit mines.

Gibbs sampling is routinely used to sample truncated Gaussian distributions. These distributions naturally occur when associating latent Gaussian fields to category fields obtained by discrete simulation methods like multipoint, sequential indicator simulation and object-based simulation. The latent Gaussians are often used in data assimilation and history matching algorithms. When the Gibbs sampling is applied on a large lattice, the computing cost can become prohibitive. The usual practice of using local neighborhoods is unsatisfying as it can diverge and it does not reproduce exactly the desired covariance. A better approach is to use Gaussian Markov Random Fields (GMRF) which enables to compute the conditional distributions at any point without having to compute and invert the full covariance matrix. As the GMRF is locally defined, it allows simultaneous updating of all points that do not share neighbors (coding sets). We propose a new simultaneous Gibbs updating strategy on coding sets that can be efficiently computed by convolution and applied with an acceptance/rejection method in the truncated case. We study empirically the speed of convergence, the effect of choice of boundary conditions, of the correlation range and of GMRF smoothness. We show that the convergence is slower in the Gaussian case on the torus than for the finite case studied in the literature. However, in the truncated Gaussian case, we show that short scale correlation is quickly restored and the conditioning categories at each lattice point imprint the long scale correlation. Hence our approach enables to realistically apply Gibbs sampling on large 2D or 3D lattice with the desired GMRF covariance. GMRF template for the 3D case is derived.Based on coding sets, exact and vastly parallel Gibbs sampling is implemented.Convergence rate of Gibbs sampler is studied.New method allows quick simulation of truncated Gaussian by acceptance.Provide latent Gaussians to represent large discrete fields for history matching.

This article describes the use of non-stationary covariance functions with compact support to estimate and simulate a random function. Based on the kernel convolution theory, the functions are derived by convolving hyperspheres in $$\mathbb{R}^n$$ followed by a Radon transform. The order of the Radon transform controls the differentiability of the covariance functions. By varying spatially the hyperspheres radius one defines non-stationary isotropic versions of the spherical, the cubic and the penta-spherical models. Closed-form expressions for the non-stationary covariances are derived for the isotropic spherical, cubic, and penta-spherical models. Simulation of the different non-stationary models is easily obtained by weighted average of independent standard Gaussian variates in both the isotropic and the anisotropic case. The non-stationary spherical covariance model is applied to estimate the overburden thickness over an area composed of two different geological domains. The results are compared to the estimation with a single stationary model and the estimation with two stationary models, one for each geological domain. It is shown that the non-stationary model enables a reduction of the mean square error and a more realistic transition between the two geological domains.