Investments in oil and gas projects are driven by critical field development decisions including well placement, which often significantly affect the projects’ economics. Due to their typically high cost and inherently insufficient data (especially in greenfields, or fields in their early stage of development), managing uncertainty is critical when optimizing a field development plan. The use of a single deterministic base case for hydrocarbon, both in place assessment and production forecasting is often misleading and leads to sub-optimal decisions. Consequently, robust field development plans require multiple geological realizations covering the range of uncertainty in reservoir properties and encompassing both multiple geological concepts and geostatistical properties distribution. Typically, an objective function such as the average net present value (NPV) or the average cumulative oil production (COP) is optimized in order to select an optimal development scenario. Nevertheless, such an assessment can be computationally prohibitive, especially when using optimization methods require hundreds, often thousands of costly simulations over a single realization, a number that significantly increases when multiple realizations are involved. This study proposes a new method for well placement optimization under uncertainty, building on map-based evolutionary optimization technique: the black hole particle swarm optimization (BHPSO). The statistical net hydrocarbon thickness (SNHCT) map is introduced to guide the BHPSO algorithm; and hence, pragmatically account for uncertainty in the process of well placement optimization. We optimize well placement on the realization corresponding to the minimum difference between its NHCT map and the SNHCT map. The SNHCT combines the average and the P90 NHCT maps; hence, assuring that the selected sweet spots for well placement are statistically the best with regard to the multiple subsurface realizations. The method is applied on the Olympus benchmark case and results are compared to two scenario reduction methods: RMfinder and k-means-k-medoids Clustering. Results show superior performance over the two methods in terms of optimality of the result and the required computational load.

Ensemble data assimilation methods, particularly iterative forms of ensemble smoother, are very useful assisted history matching techniques. One of the main challenges in the application of these methods is the excessive variance loss due to sampling errors occasioned by the limited size ensembles used in practical cases. The standard procedure to mitigate this problem is called distance-based localization. However, in this case, the data assimilation performance becomes highly dependent on the choice of the localization region. Moreover, there are several relevant examples of non-local model parameters and data types, in which cases, distance-based localization does not apply. This paper proposes a distance-free localization procedure that combines a least-squares support vector (LS-SVR) proxy with a non-isotropic tapering function proposed in the literature. We tested the proposed method in two versions of the PUNQ-S3 problem. In the first version, we tested the method to localize scalar (non-local) parameters. The results show relevant improvements, especially in terms of preserving the posterior variance of the ensemble compared to methods investigated in Lacerda et al. (J. Pet. Sci. Eng. 172,690–706, 2019). The second version considered the problem of localization of grid parameters. In this case, the proposed method outperformed a distance-based localization approach proposed in the literature.

In this work, we use the method of volume averaging to upscale the pore-scale diffusion equation on the Sierpinski carpet. Based on the isotropy condition in the fractal structure and the fact that the ratio of length scales in the Sierpinski carpet is constant, a general expression for the effective diffusion coefficient regarding the iteration of the fractal structure is suggested. Additionally, a general expression is recommended when such a ratio is not constant. The comparison between the direct numerical simulations and the results obtained from the upscaled model suggests that using a simplified expression for the effective diffusion coefficient is an attractive option when simulating large-scale fractal systems.

Water coning is one of the common issues in subsurface systems in which water flows into the production well through perforated zones. This phenomenon can cause severe problems in wellbore and surface facilities. Thus, accurate prediction of water breakthrough can help to adapt to the production mode and avoid such issues. Conducting flow simulations, as a conventional approach, can be very time demanding if one deals with large subsurface systems. Furthermore, several types of data are often collected during the life of a subsurface system each of which can help to predict the breakthrough and water coning. As such, it is very important to produce similar results using the time sequence data gathered from various geo-sensing tools. In this paper, a deep long short-term memory (LSTM) model is developed to predict the water cut and water breakthrough time for multiple production wells in a water flooding case. The dataset is generated by the Egg model with a multi-input-multi-output system. We found that the proposed model can capture the general trend of variation for the water cut time sequence data for a complex subsurface system. To evaluate the performance of our data-driven method, the results are compared with vanilla recurrent neural network (RNN), deep gated recurrent unit (GRU), and artificial neural network (ANN). The conducted comparison indicates that the proposed deep LSTM model outperforms the other three approaches when the results are compared with the numerical data.

In this paper, a novel computational framework is introduced for simulation of multiphase flow, geomechanics, and fracture propagation in porous media based on Biot’s model for poroelasticity by focusing on interactions between hydraulic and natural fractures. Since realistic porous media contain many natural fractures, it is important not only to stimulate hydraulic fractures but also to study the interaction between natural and hydraulic fractures. Here, state-of-the-art numerical modeling of natural and hydraulic fractures using a diffusive adaptive finite element phase field approach is employed. The locally mass conservative enriched Galerkin finite element methods (EG) are utilized to model two-phase flow in propagating fractures with relative permeability and capillary pressure. Geomechanics approximated by a continuous Galerkin finite element method is coupled to multiphase flow by applying an iteratively coupled scheme. Numerical examples are presented that demonstrate the effectiveness of this framework for different propagation scenarios by varying the degrees of physics. In addition, the capabilities to perform high-fidelity simulations on complex fracture networks, with randomly joined diffusive natural fractures, are illustrated.

In this work, we consider a hybrid mixed finite element method for Biot's model. The hybrid P1-RT0-P0 discretization of the displacement-pressure-Darcy's velocity system of Biot's model presented in \cite{C. Niu} is not uniformly stable with respect to the physical parameters, resulting in some issues in numerical simulations. To alleviate such problems, following \cite{V. Girault}, we stabilize the hybrid scheme with face bubble functions and show the well-posedness with respect to physical and discretization parameters, which provide optimal error estimates of the stabilized method. We introduce a perturbation of the bilinear form of the displacement which allows for the elimination of the bubble functions. Together with eliminating Darcy's velocity by hybridization, we obtain an eliminated system whose size is the same as the classical P1-RT0-P0 discretization. Based on the well-posedness of the eliminated system, we design block preconditioners that are parameter-robust. Numerical experiments are presented to confirm the theoretical results of the stabilized scheme as well as the block preconditioners.

In this work, we consider the transport of a surfactant in variably saturated porous media. The water flow is modelled by the Richards equations and it is fully coupled with the transport equation for the surfactant. Three linearization techniques are discussed: the Newton method, the modified Picard, and the L-scheme. Based on these, monolithic and splitting schemes are proposed and their convergence is analyzed. The performance of these schemes is illustrated on five numerical examples. For these examples, the number of iterations and the condition numbers of the linear systems emerging in each iteration are presented.

Predicting the petrophysical properties of rock samples using micro-CT images has gained significant attention recently. However, an accurate and an efficient numerical tool is still lacking. After investigating three numerical techniques, (i) pore network modeling (PNM), (ii) the finite volume method (FVM), and (iii) the lattice Boltzmann method (LBM), a workflow based on machine learning is established for fast and accurate prediction of permeability directly from 3D micro-CT images. We use more than 1100 samples scanned at high resolution and extract the relevant features from these samples for use in a supervised learning algorithm. The approach takes advantage of the efficient computation provided by PNM and the accuracy of the LBM to quickly and accurately estimate rock permeability. The relevant features derived from PNM and image analysis are fed into a supervised machine learning model and a deep neural network to compute the permeability in an end-to-end regression scheme. Within a supervised learning framework, machine and deep learning algorithms based on linear regression, gradient boosting, and physics-informed convolutional neural networks (CNNs) are applied to predict the petrophysical properties of porous rock from 3D micro-CT images. We have performed the sensitivity analysis on the feature importance, hyperparameters, and different learning algorithms to make a prediction. Values of R2 scores up to 88% and 91% are achieved using machine learning regression models and the deep learning approach, respectively. Remarkably, a significant gain in computation time—approximately 3 orders of magnitude—is achieved by applied machine learning compared with the LBM. Finally, the study highlights the critical role played by feature engineering in predicting petrophysical properties using deep learning.

A common method of imaging with seismic data is reverse time migration (RTM) in which recorded data are back-propagated to form an image of the subsurface. Least-squares RTM (LSRTM) extends this method to an iterative method that minimizes a data misfit term. An appropriate regularization term such as a sparsity-promoting functional (e.g., total variation (TV)) is required to stabilize the LSRTM solution. In this paper, in order to efficiently solve such regularized LSRTM via distributed optimization algorithms, we first reformulate the problem into a consensus form. The alternating direction method of multipliers (ADMM) allows us to split the problem into separate subproblems, resulting to a two-step iteration which efficiently solves the original problem. The most time-consuming step is due to LSRTM of the data associated with each source, separately, which is carried out in parallel via a set of workers. The resulting local images are sent to a master that is responsible for generating (via a proximal mapping) a unique global image that is close to the mean of local images while minimizing the regularizing function. This loop gives a general framework for the parallel computation of regularized LSRTM and allows different regularizing functions to be employed by using appropriate proximal operators. We demonstrate the performance of the proposed algorithm with a set of numerical examples. The results show that TV methods can generate more accurate images compared with the l1-norm sparse LSRTM, l2-norm LSRTM, and traditional RTM. Furthermore, the new algorithm shows improvements over the traditional formulations due to the consensus parallelization.