The purpose of this page is to review the mathematics of fluid flow. We limit our review to essential aspects of partial differential equations, vector analysis, numerical methods, matrices, and linear algebra. These topics are discussed in the context of two fluid flow applications: analysis of the convection/dispersion equation and diagonalization of the permeability tensor. For more details about the mathematics presented here, consult the literature. Partial differential equations (PDEs) are frequently encountered in petroleum engineering.
Crandall, Dustin (National Energy Technology Laboratory) | Gill, Magdalena (National Energy Technology Laboratory, LRST) | Moore, Johnathan (National Energy Technology Laboratory, LRST) | Brown, Sarah (West Virginia Geological and Economic Survey) | Mackey, Paige (National Energy Technology Laboratory, ORISE)
The behavior of fractured low-permeability rock in many subsurface formations is critical for unconventional resource extraction. Understanding how flow through individual fractures changes during shearing, and what influence heterogeneity of the rock has on shearing behavior, was the focus of our laboratory study. Computed tomography (CT) scanning of fractured rocks undergoing shear was coupled with numerical simulations of fluid flow through these fractures. We sheared multiple cores from the Marcellus and Eau Claire shales in a closed system with confining pressures of greater than 1000 psi. Samples were manually sheared in a step wise fashion. After each shearing event we assessed the bulk hydrodynamic response by measuring permeability through the core and performed a high-resolution CT scan to understand how the principal and secondary fractures were changing in the core volume. The mineralogy of each sample was examined via x-ray fluorescence.
A range of interdependent characteristics influence fracture network evolution and sample cohesion: mineralogy, lithological heterogeneity, principal fracture morphology, fracture asperities, and shearing direction in relation to bedding. We found that samples sheared parallel to bedding were less likely to develop extensive networks of secondary fractures, with secondary fracture growth contingent on the presence of large asperities. Fracture permeability tended to increase with continued shear and secondary fracture development, but a high variance existed between samples. In some instances, permeabilities decreased in response to shear-initiated aperture reduction due to fracture mating. Gouge formation is another factor contributing to the transmissivity decreases, particularly in shale-dominated fracture regions. The ability to study this complex behavior in a controlled fashion using CT scanning enables a view into processes that impact production in many unconventional formations. Findings show that small scale features and details can play a significant role in fracture behavior and should be accounted for.
Shale properties vary significantly and understanding how fractures evolve due to geomechanical stressing can improve our understanding of how to effectively stimulate a variety of formations.While hydraulic fracturing is a large-scale activity, the microfabric and heterogeneity of shale can control fracture evolution and flow properties. Upscaling the impact of microfabric and heterogeneity is poorly captured in most modeling and planning efforts; this disconnect between small scale features and large-scale operations is understandable. It is difficult to measure changes in fractures directly, difficult to implement upscaled equations of value, and difficult to know if studied laboratory/outcrop samples are representative of activities in the subsurface. This study describes the observed behavior of two distinctly different shales under controlled geomechanical stressing to examine what impact small features have on fracture evolution. By examining two shales with distinctly different structure and composition our goal is to understand when inclusions of micro-features in upscaling is critical to understanding system dynamics.
Objectives - Image-Based Rock Physics (IBRP) simulation of petrophysical properties based on sub-micron to micron-scale images of very fine-grain rocks is constrained by the resolution and range of various imaging techniques used. Unlike some conventional sandstone and carbonate reservoir rock where a single Micro-Computed X-ray Tomography (μCT) volume images nearly all of the significant pores and pore throats, many low-permeability rock types contain phase regions with micro-pores and pore throats, including intergranular microcrack pores, that are not accurately resolved at the required μCT scale needed for a representative elementary volume (REV) for the whole rock. Properties for these regions are obtained at a finer-scale or using a different measurement method and these properties then assigned to the phase regions at the larger REV scale. This study explores the methodology involved in obtaining and assigning microcrack properties in μCT rock images and demonstrates a workflow to handle uncertainty in the location and properties of microcracks using two representative low-permeability sandstones.
Methods/Procedures/Process - The workflow combines μCT images of a mini-plug sample (~50mm3), which represents the rock REV, with Focused Ion Beam - Scanning Electron Microscopy (FIB-SEM) images (~200μm3) of regions of various types of observed microporosity (including intergranular microcrack pores) which occur within the REV sample. Different representative types of microporosity regions were imaged and properties calculated from the higher-resolution FIB-SEM image volumes. For some fraction of μCT microporosity regions, such as micro-fractures, their locations in the REV μCT sample was known but the micro-fracture properties were not known. A sub-resolution micro-fracture model was numerically constructed, honoring the mineral facies morphology and microporosity types assigned based on their respective distributions as observed in high resolution SEM images. Resultant porosity, capillary pressure and flow properties on the larger REV volume were cross-validated with independent core analysis measurements.
Results/Observations/Conclusions - This study illustrates a workflow for assigning properties, obtained at finer scales or using other measurement methods, to regions in the REV at larger scale but lower resolution. The resulting rock model produces the same porosity, permeability, and capillary pressure as core analysis measurements, and has the potential to predict relative permeability.
Applications/Significance/Novelty - It is expected that the majority of low-permeability rocks require an upscaling methodology similar to that developed in this study for IBRP computations and integration with core analysis. Using this methodology IBRP offers deeper understanding of building blocks of the upscaled-properties measured by core analysis. IBRP also offers the ability to measure/compute relative permeabilities that are nearly physically impossible to measure on core and the ability to construct digital rocks that allow evaluation of complete suites of rocks and their properties.
Zhang, Na (Division of Sustainable Development, College of Science and Engineering, Hamad Bin Khalifa University) | Abushaikha, Ahmad Sami (Division of Sustainable Development, College of Science and Engineering, Hamad Bin Khalifa University)
Modelling fluid flows in fractured reservoirs is crucial to many recent engineering and applied science research. Various numerical methods have been applied, including finite element methods, finite volume methods. These approaches have inherent limitations in accuracy and application. Considering these limitations, in this paper, we present a novel mimetic finite difference (MFD) framework to simulate two phase flow accurately in fracture reservoirs.
A novel MFD method is proposed for simulating multiphase flow through fractured reservoirs by taking advantage of unstructured mesh. Our approach combines MFD and finite volume (FV) methods. Darcy's equation is discreted by MFD method, while the FV method is used to approximate the saturation equation. The resulting system of equations is then imposed with suitable physical coupling conditions along the matrix/ fracture interfaces. This coupling conditions at the interfaces between matrix and fracture flow involve only the centroid pressure of fractures, which brings some simplification in analysis. The proposed approach is applicable for three dimensional systems. Moreover, it is applicable in arbitrary unstructured gridcells with full-tensor permeabilities. Some examples are implemented to show the performance of MFD method. The results showed a big potential of our method to simulate the flow problems with high accuracy and application.
Ahmadinia, Masoud (Centre for Fluid and Complex Systems, Coventry University) | Shariatipour, Seyed (Centre for Fluid and Complex Systems, Coventry University) | Andersen, Odd (SINTEF Digital, Mathematics and Cybernetics) | Sadri, Mahdi (Centre for Fluid and Complex Systems, Coventry University)
To improve the reservoir simulation model, uncertain parameters such as porosity and permeability in the reservoir rock strata need to be adjusted to match the simulated production data with the actual production data. This process is known as History Matching (HM). In geological CO2 storage that is being promoted for use in depleted hydrocarbon reservoirs and saline aquifers, CO2 tends to migrate upwards and accumulate as a separate plume in the zone immediately beneath the reservoir caprock. Thus caprock morphology is of considerable importance with respect to storage safety and migration prediction for the purpose of long-term CO2 storage. Moreover, small scale caprock irregularities, which are not captured by seismic surveys, could be one of the sources of errors while matching the observed CO2 plume migration and the numerical modelling results (e.g. Sleipner). Thus here we study the impact of uncertainties in slope and rugosity (small scale caprock irregularities not captured by seismic surveys) on plume migration, using a history-matching process. We defined 10 cases with different initial guesses to reproduce the caprock properties representing an observed plume shape. The results showed a reasonable match between the horizontal plume shape in calibrated and observed models with an average error of 2.95 percentages
Algebraic Multigrid (AMG) methods have proven to be efficient when numerically solving elliptic Partial Differential Equations (PDE). In reservoir simulation, AMG is used together with the Constrained Pressure Residual (CPR) method to solve a partially decoupled pressure system. Recently, effort has been focused on improving the robustness of the AMG-CPR solver. This paper presents the performance of different AMG-CPR strategies for massive reservoir models. In addition a solver selection analysis is conducted proving that dynamic selection of solvers has potential of increasing the overall efficiency and robustness of the simulation.
Numerous decoupling/preconditioning algorithms exist and have been shown to influence the pressure matrix properties, some resulting in matrices more suitable to the characteristics favorable to AMG. Several decoupling/preconditioning strategies are investigated such as Alternate Block Factorization (ABF), Quasi-IMPES (QI), and Dynamic Rowsum (DRS). The extracted pressure matrix could be suitable or unsuitable for AMG, depending on the matrix row sum, the diagonal signs, and the signs of the off-diagonal values.
The advantage of using AMG as a preconditioner is demonstrated by running the SPE10 case. The recommended AMG settings that result in the optimal performance for SPE10 are shared. A speedup is seen of up to 4X when using AMG with optimal settings versus the default solver in the in-house reservoir simulator with the improvement range depending on the number of processors used. SPE10 is a highly heterogeneous model resulting in matrices favorable for AMG, i.e. pressure decoupling produces positive definite pressure matrices which is not necessarily representative of industry models. A comparison is then made with a selection of models with a wide range of characteristics and finally an examination of the convergence behavior of key industry cases with different decoupling strategies is presented. The overall convergence behavior of the pressure and full system are shown and the top decoupling algorithms for the particular models are discussed. Finally, the applicability and performance gain of selectively using AMG during a run is demonstrated.
Recent developments have been made in regard to AMG methods; but their applicability in a wide range of massive real cases is yet to be explored. In this work, different decoupling methods are tested, the AMG behavior on real field massive models is analyzed, the scalability is investigated, and AMG is selectively activated during a simulation run shedding light on the potential of future work entailing the use of Artificial Intelligence (AI) to dynamically select the optimal solver choice.
Abd, Abdul Salam (Division of Sustainable Development College of Science and Engineering Hamad Bin Khalifa University) | Abushaikha, Ahmad (Division of Sustainable Development College of Science and Engineering Hamad Bin Khalifa University)
Reactive transport is an area of growing interest to the petroleum industry due to the need to develop efficient simulators for oil production in mature and fractured fields. The subsurface flow processes are highly dependent on the chemical reaction between fluid and rocks, and between fluids themselves. This means that fluid flow in a reservoir should be characterized by mass change and chemical reactions, and thus researchers are trying to account for both phenomena through their corresponding governing equations. In this work, we discuss the early efforts to couple mass transfer and chemical reaction equations through presenting key studies in the field of modelling reactive transport and suggest a new discretization scheme for geochemical reactions. Our main objective is to improve on the formulation and simulation of reactive transport problems and explain how to implement Mimetic Finite-Difference (MFD) discretization scheme in reservoir simulators. The paper outlines the steps to discretize the governing equations of reactive transport, and construct and solve the resultant Jacobian system. The purpose of this work is to provide a clear presentation of the main mathematical and theoretical concepts of reactive transport, and how they can be applied in reservoir simulation framework. Such framework can be utilized in the future to develop a state-of-art reservoir simulator that employs Mimetic Finite Difference schemes in unstructured grids and full tensor permeability structures to solve for fluid flow in porous media, while accounting for the geochemical reactions at the subsurface.
Da Silva Moreira, Paulo Henrique (LRAP / Universidade Federal do Rio de Janeiro) | Gomes da Silveira, Thaís Márcia (LRAP / Universidade Federal do Rio de Janeiro) | Drexler, Santiago (LRAP / Universidade Federal do Rio de Janeiro) | Couto, Paulo (LRAP / Universidade Federal do Rio de Janeiro)
The reliable prediction of reservoir performance requires the cost effective implementation of oil recovery systems, and it is necessary to simulate the fluid flow processes in the reservoir and to measure the rock and fluid properties that determine reservoir behaviour. However, a good prediction relies on accurate values of reservoir physical properties. Carbonates rocks in Brazilian Pre-salt are known for their heterogeneity. Characterizing their physical represents a great challenge and the combination of experimental and computational techniques lead to a more comprehensive understanding of the reservoir behavior.
In the present work, the relative permeability curves of a carbonate core sample with respect to oil and water are calculated by matching the data obtained in a labscale unsteady-state core flood experiment carried out at high pressure high temperature characteristics of Brazilian Pre-salt reservoirs. Corey-type equations were used to model the relative permeability due to its simplicity and having fewer parameters involved. The Monte Carlo Markov chain (MCMC) method was used as optimization tool, taking the fluid production and pressure drop measurements collected during the core flood experiment as input data. An alalysis of the sensitivity cofficients was carried out in order to deal with eventual linear dependences among the terms to be estimated. The Markov chain was generated and its convergence observed. The posterior distributions of the constant terms in the Corey equations were calculated and their mean values applied in order to calculate the relative permeability curves for the oil and water phases. The range of water saturation in which the relative permeability curves describe the core conditions after the breakthrough time, due to the occurrence of capillary end-effects, was calculated. The history match of fluid production and pressure drop was carried out, showing a good fit between the pressure curves. A gap was observed between the production curves due to the fact that the experimental measurements accounted the cumulative volume of oil and water, while the theoretical curve accounted the oil volume only.
Immiscible water-alternating-gas (iWAG) flooding is often considered as a tertiary recovery technique in waterflooded or about-to-be waterflooded reservoirs to increase oil recovery due to better mobility control and potentially favorable hysteretic changes to phase relative permeabilities. In such cases, typically, reservoir simulation models already exist and have been calibrated, often modifying saturation functions during the history matching stage. However, to utilize such models in forecasting iWAG performance, additional parameters may be required. These can be acquired by simulation of WAG coreflood experiments. While in many published cases, the parameter values obtained from matching experimental results are used without modification, this may not be advisable since the parameters are only valid at the core scale at which they were obtained. This paper discusses the challenge of systematically upscaling WAG parameters obtained at core scale to an existing full field model.
In this work, we use a multi-stage upscaling process from core scale to full field scale. The first stage uses a core scale model to match ‘representative’ core flood experiments and obtain WAG parameters. The second uses a well-to-well high-resolution 1D section of the full field model populated using gridblocks of core size to generate ‘reference’ WAG performance using the unaltered WAG parameters obtained from core. The third stage uses a similar 1D model but populated using gridblocks at full field model resolution to match the results from the reference model while adjusting the WAG parameters as little as possible. Finally, a model using the full field model resolution as well as the full field relative permeability functions which, it is assumed, have been tuned to match the history and account for dispersion is used to match the reference model results and obtain final upscaled WAG parameters.
The upscaled WAG parameters obtained at the end of this multi-stage process can be used at the field scale. This process allows clear quantification of the uncertainty associated with the upscaling process. Simulations at the third stage showed that once the full field to core scale grid size ratio exceeded a certain point (2500:1), there was a marked increase in the difference between upscaled and reference model results. It was found that if WAG parameters were changed in the full field model resolution model in order to match recovery results in the reference model, Land's parameter could change by up to 10% and relative permeability reduction factor could increase by up to 30% although it is expected that this will vary from case to case. It is therefore recommended to identify and use full field model resolutions to as close to the threshold as possible. The practice of using the core scale iWAG parameters in the full field model directly could under-estimate actual recovery, and overestimate injectivity. When considering the WAG mechanism alone, the value of the recovery underestimate increasing with pore volumes injected and, in our case, by up to 7% after injecting 1 pore volume of fluid.
This multi-stage simulation approach helps identify the adjustments required and uncertainties associated with simulating iWAG flooding in reservoir models. This approach utilizes options widely present in commercially available finite difference simulators, addresses the challenge of utilizing existing pseudo functions and provides a practical methodology through which iWAG performance forecasting can be improved.
Sun, Zheng (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing)) | Shi, Juntai (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing)) | Wu, Keliu (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing)) | Zhang, Tao (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing)) | Feng, Dong (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing)) | Li, Xiangfang (MOE Key Laboratory of Petroleum Engineering and State Key Laboratory of Petroleum Resources and Engineering, China University of Petroleum (Beijing))
Low-permeability coalbed-methane (CBM) reservoirs possess unique pressure-propagation behavior, which can be classified further as the expansion characteristics of the drainage area and the desorption area [i.e., a formation in which the pressure is lower than the initial formation pressure and critical-desorption pressure (CDP), respectively]. Inevitably, several fluid-flow mechanisms will coexist in realistic coal seams at a certain production time, which is closely related to dynamic pressure and saturation distribution. To the best of our knowledge, a production-prediction model for CBM wells considering pressure-propagation behavior is still lacking. The objective of this work is to perform extensive investigations into the effect of pressure-propagation behavior on the gas-production performance of CBM wells. First, the pressure-squared approach is used to describe the pressure profile in the desorption area, which has been clarified as an effective-approximation method. Also, the pressure/saturation relationship that was developed in our previous research is used; therefore, saturation distribution can be obtained. Second, an efficient iteration algorithm is established to predict gas-production performance by combining a new gas-phase-productivity equation and a material-balance equation. Finally, using the proposed prediction model, we shed light on the optimization method for production strategy regarding the entire production life of CBM wells. Results show that the decrease rate of bottomhole pressure (BHP) should be slow at the water single-phase-flow stage, fast at the early gas/water two-phase-flow stage, and slow at the late gas/water two-phase-flow stage, which is referred to as the slow/fast/slow (SFS) control method. Remarkably, in the SFS control method, the decrease rate of the BHP at each period can be quantified on the basis of the proposed prediction model. To examine the applicability of the proposed SFS method, it is applied to an actual CBM well in Hancheng Field, China, and it enhances the cumulative gas production by a factor of approximately 1.65.