Thiele, Christopher (Rice University) | Araya-Polo, Mauricio (Shell International Exploration & Production, Inc.) | Alpak, Faruk Omer (Shell International Exploration & Production, Inc.) | Riviere, Beatrice (Rice University)
Direct numerical simulation of multiphase pore-scale flow is a computationally demanding task with strong requirements on time-to-solution for the prediction of relative permeabilities. In this paper, we describe the hybrid-parallel implementation of a two-phase two-component incompressible flow simulator using MPI, OpenMP, and general-purpose graphics processing units (GPUs), and we analyze its computational performance. In particular, we evaluate the parallel performance of GPU-based iterative linear solvers for this application, and we compare them to CPUbased implementations of the same solver algorithms. Simulations on real-life Berea sandstone micro-CT images are used to assess the strong scalability and computational performance of the different solver implementations and their effect on time-to-solution. Additionally, we use a Poisson problem to further characterize achievable strong and weak scalability of the GPU-based solvers in reproducible experiments. Our experiments show that GPU-based iterative solvers can greatly reduce time-to-solution in complex pore-scale simulations. On the other hand, strong scalability is currently limited by the unbalanced computing capacities of the host and the GPUs. The experiments with the Poisson problem indicate that GPU-based iterative solvers are efficient when weak scalability is desired. Our findings show that proper utilization of GPUs can help to make our two-phase pore-scale flow simulation computationally feasible in existing workflows.
Araujo, Mariela (Shell International Exploration and Production Inc.) | Chen, Chaohui (Shell International Exploration and Production Inc.) | Gao, Guohua (Shell International Exploration and Production Inc.) | Jennings, Jim (Shell International Exploration and Production Inc.) | Ramirez, Benjamin (Shell International Exploration and Production Inc.) | Xu, Zhihua (ExxonMobil) | Yeh, Tzu-hao (Shell International Exploration and Production Inc.) | Alpak, Faruk Omer (Shell International Exploration and Production Inc.) | Gelderblom, Paul (Shell International Exploration and Production Inc.)
Increased access to computational resources has allowed reservoir engineers to include assisted history matching (AHM) and uncertainty quantification (UQ) techniques as standard steps of reservoir management workflows. Several advanced methods have become available and are being used in routine activities without a proper understanding of their performance and quality. This paper provides recommendations on the efficiency and quality of different methods for applications to production forecasting, supporting the reservoir-management decision-making process.
Results from five advanced methods and two traditional methods were benchmarked in the study. The advanced methods include a nested sampling method MultiNest, the integrated global search Distributed Gauss-Newton (DGN) optimizer with Randomized Maximum Likelihood (RML), the integrated local search DGN optimizer with a Gaussian Mixture Model (GMM), and two advanced Bayesian inference-based methods from commercial simulation packages. Two traditional methods were also included for some test problems: the Markov-Chain Monte Carlo method (MCMC) is known to produce accurate results although it is too expensive for most practical problems, and a DoE-proxy based method widely used and available in some form in most commercial simulation packages.
The methods were tested on three different cases of increasing complexity: a 1D simple model based on an analytical function with one uncertain parameter, a simple injector-producer well pair in the SPE01 model with eight uncertain parameters, and an unconventional reservoir model with one well and 24 uncertain parameters. A collection of benchmark metrics was considered to compare the results, but the most useful included the total number of simulation runs, sample size, objective function distributions, cumulative oil production forecast distributions, and marginal posterior parameter distributions.
MultiNest and MCMC were found to produce the most accurate results, but MCMC is too costly for practical problems. MultiNest is also costly, but it is much more efficient than MCMC and it may be affordable for some practical applications. The proxy-based method is the lowest-cost solution. However, its accuracy is unacceptably poor.
DGN-RML and DGN-GMM seem to have the best compromise between accuracy and efficiency, and the best of these two is DGN-GMM. These two methods may produce some poor-quality samples that should be rejected for the final uncertainty quantification.
The results from the benchmark study are somewhat surprising and provide awareness to the reservoir engineering community on the quality and efficiency of the advanced and most traditional methods used for AHM and UQ. Our recommendation is to use DGN-GMM instead of the traditional proxy-based methods for most practical problems, and to consider using the more expensive MultiNest when the cost of running the reservoir models is moderate and high-quality solutions are desired.
Guohua Gao and Jeroen C. Vink, Shell Global Solutions US; and Faruk O. Alpak, Shell International Exploration and Production Summary The in-situ upgrading process (IUP) is a thermal-recovery technique that relies on a pattern-based development process, a complicated physical process that involves thermal and mass transfer in porous media, which renders full field-scale reservoir simulations impractical. Although it is feasible to quantify the impact of subsurface uncertainties on recovery for small-scale sector models with experimental design (ED), it is still a very challenging problem to quantify their impact on field-scale quantities. Straightforward upscaling to field scale does not work because such conventional superposition-based methods do not capture the effects of spatial variability in rock and fluid properties and the time delay in sequential pattern development. In this paper, we show that, under certain mild assumptions, an analytical superposition formulation can be developed that propagates the uncertainties of production forecasts and economic evaluations generated from a sector model to full field-scale quantities. One can simplify this formulation further so that the variance of a field-scale quantity is analytically expressed as the variance of the same single-pattern quantity multiplied by a (computable) scaleup factor. This makes it possible to implement a practical uncertainty quantification work flow in which single-pattern results are upscaled to accurate full field results with reliable uncertainty ranges, without the need for full field-scale simulations. We apply the proposed novel superposition and uncertaintypropagation method to a multipattern IUP development, and demonstrate that this work flow produces reliable results for field-scale production and economics as well as realistic uncertainty ranges. Moreover, these results indicate that the scaleup factor for singlepattern results can accurately capture the impact of spatial correlations of subsurface uncertainties, the size of the field-scale model, the time-delay in pattern development, and the discount rate.
Representing the complete spectrum of fine-scale stratigraphic details in full-field dynamic models of geologically complex clastic reservoirs is beyond the limits of existing computational capabilities. A quasiglobal multiphase upscaling method--the regional-scale multiphase upscaling (RMU) method--is developed, in which the dynamic effects of subgrid-scale (typically subseismic) nonlocal stratigraphic reservoir elements (e.g., channels, lobes, sand bars, and shale drapes) are captured by means of pseudofunctions for flow simulation. Unlike conventional dynamic multiphase upscaling methods, the RMU method does not require fine-resolution reservoir-scale simulations. Rather, it relies on intermediate-scale sector-model simulations for pseudoization. The intermediate scale, also referred to as the regional scale, is defined as the spatial scale at which the global multiphase flow effects of nonlocal stratigraphic elements can be approximated by fine-resolution flow simulations with reasonable accuracy. During the pseudoization process, dynamic multiphase flow responses of coarse regional-scale sector models are calibrated against those stemming from their corresponding fine-resolution parent models. Each regional-scale sector model is simulated only once at the fine geologic resolution. The process involves automatic determination and subsequent modification of the parameters that describe rock relative permeability and capillary pressure functions. Coarse regional-scale models are simulated a few times until a reasonable match between their coarse- and fine-resolution dynamic responses can be attained. The parameter-estimation step of the pseudoization process is performed by use of a very efficient constrained nonlinear optimization algorithm. The RMU method is evaluated in two proof-of-concept numerical examples involving a plethora of turbidite stratigraphic architectures. The method yields simulation results that are always more accurate than conventionally upscaled coarse-resolution model predictions. Incorporating geologically based pseudofunctions into otherwise simple coarse-resolution full-field reservoir models reduces the simulation cycle time significantly and improves the accuracy of production forecasts. The RMU method typically delivers two to three orders of magnitude in speed up of flow simulations.
Chen, Tianhong (Shell International Exploration and Production Inc.) | Noirot, Jean-Christophe (Shell Nigeria E&P Co. Ltd.) | Khandelwal, Arpit (Shell India Markets Private Ltd) | Xue, Guangri (Shell Oil Co.) | Barton, Mark (Shell Intl E&P) | Alpak, Faruk Omer (Shell Intl. E&P Co.)
Well test analysis in turbidite reservoirs is complicated by the intricate stratigraphy prevailing in this depositional environment. Because of this complexity, important reservoir architectural parameters driving flow behavior (e.g., shale drape coverage, object dimensions) cannot be estimated using simple analytical reservoir models employed in conventional well test analysis techniques. Alternatively, simulation-based well test analysis offers the advantage of being able to capture stratigraphic complexity. However, it requires a very large number of models and simulations to identify multiple solutions to such a highly non-unique inversion problem. In this work, we have developed a novel well test analysis workflow by constructing a large library of build-up type curves derived by appropriately scaling a comprehensive set of reference drawdown simulations. This set is used to rapidly identify a variety of stratigraphic scenarios matching a given well test. Key stratigraphic parameters are then estimated through statistical analysis of the results. The proposed well test analysis technique has been applied to synthetic and field examples. For the tested cases, stratigraphic interpretations derived from well tests are found to be consistent with those obtained from other data sources.
Jin, Long (Shell) | Gao, Guohua (Shell) | Vink, Jeroen C. (Shell Intl. E&P Co.) | Chen, Chaohui (Shell International EP) | Weber, Daniel (Shell Intl. E&P Co.) | Alpak, Faruk Omer (Shell Intl. E&P Co.) | van den Hoek, Paul (Shell) | Pirmez, Carlos (Shell Intl. E&P Co.)
Quantitative integration of 4D seismic data with production data into reservoir models is a challenging task. One important issue is how to properly quantify the uncertainty, or the posterior probability distribution (PPD). The Very Fast Simulated Annealing (VFSA) is a stochastic searching method, whereas the Simultaneous Perturbation and Multivariate Interpolation (SPMI) is a model-based local searching method. The stochastic features of the VFSA provide the feasibility of identifying possible multiple peaks of a PPD, but it converges very slowly. On the other hand, the model-based SPMI method has the advantages of effectively utilizing the smooth features of an objective function, and thus can converge to local optimum very quickly. More importantly, the Hessian of the objective function, or the covariance matrix of the PPD, can be estimated by the SPMI method with satisfactory accuracy. However, it is very difficult to identify multiple optima by applying the SPMI method alone. In this paper, we propose an efficient joint inversion workflow by appropriately integrating the two derivative
free optimization (DFO) methods. The complementary features of the two methods can further improve both applicability and efficiency of this joint inversion workflow. We tested the workflow with a 3D synthetic model and a real field case. Our results show that the integrated method is efficient and can deliver good results for jointly assimilating 4D seismic and production data.
Chen, Chaohui (Shell International EP) | Jin, Long (Shell) | Gao, Guohua (Shell) | Weber, Daniel (Shell Intl. E&P Co.) | Vink, Jeroen C. (Shell Intl. E&P Co.) | Hohl, Detlef (Shell Intl. E&P BV) | Alpak, Faruk Omer (Shell Intl. E&P Co.) | Pirmez, Carlos (Shell Exploration and Production Company)
Gradient-based optimization algorithms can be very efficient in history matching problems. Since many commercial reservoir simulators do not have an adjoint formulation built in, exploring capability and applicability of derivative-free optimization (DFO) algorithms is crucial. DFO algorithms treat the simulator as a black box and generate new searching points using objective function values only. DFO algorithms usually require more function evaluations, but this obstacle can be overcome by exploiting parallel computing.
This paper tests three DFO algorithms, Very Fast Simulated Annealing (VFSA), Simultaneous Perturbation and Multivariate Interpolation (SPMI) and Quadratic Interpolation Model-based (QIM) algorithm. Both SPMI and QIM are model-based methods. The objective function is approximated by a quadratic model interpolating points evaluated in previous iterations, and new search points are obtained by minimizing the quadratic model within a trust region. VFSA is a stochastic search method. These algorithms were tested with two synthetic cases (IC fault model and Brugge model) and one deepwater field case. Principal Component Analysis is applied to the Brugge case to parameterize the reservoir model vector to less than 40 parameters.
We obtained good matches with all three derivative-free methods. In terms of number of iterations used for converging and the final converged value of the objective function, SPMI outperforms the others. Since SPMI generates a large number of perturbation and search points simultaneously in one iteration, it requires more computer resources. QIM does not generate as many interpolation points as SPMI, and it converges more slowly in terms of time. VFSA is a sequential method and usually requires hundreds of iterations to converge.
Jin, Long (Shell) | Alpak, Faruk Omer (Shell Intl E&P Co) | van den Hoek, Paul (Shell) | Pirmez, Carlos (Shell Intl E&P Co) | Fehintola, Tope (Shell U.K. Limited) | Tendo, Fidelis (Shell Nigeria E&P) | Olaniyan, Elozino Enite
A robust and efficient simulation technique is developed based on the extension of the mimetic finite volume method to multiscale hierarchical hexahedral (corner-point) grids via use of the multiscale mixed finite element method. The implementation of the mimetic subgrid discretization method is compact and generic for a large class of grids, and thereby, suitable for discretizations of reservoir models with complex geologic architecture. Flow equations are solved on a coarse grid where basis functions with subgrid resolution accurately account for subscale variations from an underlying fine-scale geomodel. The method relies on the construction of approximate velocity spaces that are adaptive to the local properties of the differential operator. A variant of the method for computing velocity basis functions is developed that utilizes an adaptive local-global algorithm to compute multiscale velocity basis functions by capturing the principal characteristics of global flow. Both local and local-global methods generate subgrid-scale velocity fields that reproduce the impact of fine-scale stratigraphic architecture. By using multiscale basis functions to discretize the flow equations on a coarse grid, one can retain the efficiency of an upscaling method, while at the same time produce detailed and conservative velocity fields on the underlying fine grid.
The accuracy and efficacy of the multiscale method is compared to that of fine-scale models and of coarse-scale models with no subgrid treatment for several two-phase flow scenarios. Numerical experiments involving two-phase incompressible flow and transport phenomena are carried out on high-resolution corner-point grids that explicitly represent example stratigraphic architectures found in real-life shallow marine and turbidite reservoirs. The multiscale method is several times faster than the direct solution of the fine-scale problem and yields more accurate solutions than coarse-scale modeling techniques that resort to explicit effective properties. The accuracy of the multiscale simulation method with adaptive local-global velocity basis functions are compared to that of the local velocity basis functions. The multiscale simulation results are consistently more accurate when the local-global method is employed for computing the velocity basis functions.
Multiscale problems are common in nature, especially in Earth's subsurface, a principal environment-of-focus for the hydrocarbon exploration and production industry. Multiple spatial and temporal scales exist for physical phenomena in the subsurface, from pore to core, core to well log, well log to geological model, and geological model to reservoir simulation (dynamic) model scales. The objective of our work is to develop and test the accuracy and efficiency of a novel multiscale method to bridge the gap in gridblock sizes between typical geological models (fine scale ? centimeters to decimeters in the vertical direction and meters to tens of meters in horizontal directions) and dynamic models (coarse scale ? meters to tens of meters in the vertical direction and tens of meters to hundreds of meters in horizontal directions). High degree of heterogeneity in petrophysical properties (e.g., permeability and porosity), elements of structural framework (e.g., faults, joints, and deformation bands), and stratigraphic architectural characteristics (e.g., channels, lobes, clinoforms, shale/mud drapes) of subsurface reservoirs together or individually can give rise to complex and typically large (in terms of the number of gridblocks) geological, and hence, dynamic models.
Economic constraints impose stringent limits on the number of wells that can be drilled in deepwater developments. Thus, optimal placement and operation of wells have a major impact on the project rewards. Well-placement in deepwater developments is a challenging optimization problem. Manual approaches to its solution can be cumbersome even with good use of engineering judgment: (a) There often exist many combinations of well locations subject to investigation; (b) There is need to optimize operational constraints for every well-placement scenario; (c) The optimization process has to be repeated for a variety of geologic realizations; (d) Presence of complex sub-seismic geologic architecture may render workflows that solely rely on seismic data obsolete. We developed an adjoint-based optimization algorithm that rapidly identifies alternative optimal well-placement scenarios for a given geologic realization. Adjoint-based gradients approximate the sensitivities of a suitable objective function with respect to well locations. These sensitivities guide the iterative search with improving directions that progressively maximize the ultimate recovery. The main advantage of the adjoint method is that it provides sensitivities with only one forward (reservoir) and one backward (adjoint) simulation, rendering computational costs affordable for field applications. The adjoint method in our implementation operates in local search mode.
Our algorithm is applied to two different geologically complex channelized turbidite reservoirs to identify alternative optimal locations for a new injector-producer pair given existing operational and well constraints. Ensuing results demonstrate that the adjoint-based sensitivities can be effectively used to find optimal well locations. In the field simulation models, the optimization derived well locations lead to a higher ultimate oil recovery compared to a manual optimization approach. The algorithm yields multiple alternative well-placement scenarios within a shorter amount of time compared to a single optimal location found through manual optimization.