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.
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.
Traditionally, preconditioners are used to damp slowly varying error modes in the linear solver stage. State-of-the-art multilevel preconditioners use a sequence of aggressive restriction, coarse-grid correction and prolongation operators to handle low-frequency modes on the coarse grid. High-frequency errors are then resolved by employing a smoother on fine grid. In this paper, the algebraic smoothing aggregation two level preconditioner is implemented to solve different coupled problems.
The proposed method generalizes the existing MsRSB and smoothing aggregation AMG methods. This method does not require any coarse partitioning and, hence, can be applied to general unstructured topology of the fine scale. Inspired by smoothing aggregation algebraic multigrid solver, the algebraic smoothing aggregation preconditioner constructs basis functions which allow mapping of some high-frequency modes from fine scale to low-frequency modes on the coarse scale. These basis functions are also used to reconstruct unknown primary variables at the fine scale using their approximations at the coarse level.
The proposed preconditioner has been adopted to challenging multiphysical problems, including fully coupled simulation of filtration and geomechanics processes including non-isothermal fluid flow problems. The preconditioner provides a reasonably good approximation to the coupled physical processes and speeds up the convergence. Compared to traditional ILU0+GMRES linear solvers, our preconditioner with GMRES solver reduces the number of iterations by about 3 times. In addition, the proposed method obeys a good theoretical scalability essential for parallel simulations.
Geological/flow properties computed from a heterogeneous rock sample represent the properties of the entire medium if the size of the sample is larger than the size of the representative elementary volume (REV) of the rock. Among these properties, relative permeability is believed to have the most influence on multiphase-flow behavior. Therefore, much care must be taken to ensure that relative permeability has been computed beyond the REV size. However, in the literature, experimental and numerical evaluations of relative permeability are mostly performed without such care, and the saturation functions are calculated at different flow rates while the length of the sample remains the same.
In this work, a series of numerical simulations of waterflooding are performed in a model made of alternating parallel layers of two rock types. Upscaled relative permeabilities are computed on samples at different lengths, changing from millimeter to kilometer scale. Flow behavior for different samples is analyzed according to capillary, viscous, and gravity effects, which are monitored during the simulation.
Although the REV size for some properties such as porosity seems to be the smallest replica block (unit block), our results show that the REV for relative permeability follows a power-law trend and changes with flow rate and scale. More importantly, for the first time, it is shown that the REV size is not necessarily identical for different phases. It is also demonstrated that especially at small scales, breakthrough can occur at layers where it was not expected. This study identifies the range of the length scale for which either experimental or numerical determination of relative permeability is valid.
ABSTRACT: Understanding the poro-mechanical properties of tight rocks such as shale sediments is critical and relevant in many disciplines including petrophysics, hydrology, and subsurface engineering. Due to the presence of various sources of heterogeneities, spanning across multiple length scale, these formations are suitable for multiscale modelling framework. As one of the source of heterogeneity, microcrack growth serve as an important mechanism in the inelastic and fracture behavior at the macroscale. To this end, a probabilistic four-level multiscale framework is utilized for computation of homogenized and damaged behavior at macro scale propagating uncertainty at different length scale. A micromechanics based damage model is presented within the framework of Linear Elastic Fracture Mechanics (LEFM) in which the damage evolution is represented by evolution of global crack density parameter for isotropic distribution of microcracks in transversely isotropic medium. A Monte Carlo simulation is performed to quantify the uncertainty in the model response stemming from uncertainties associated with model parameters at different length scale.
Organic-rich shale is a porous multi-scale material, consisting of clay, organic matter and silt inclusion. These formations serve as source of oil and gas, as well as the permeability barriers in a reservoir. The inherent complex multiscale and heterogeneous structure of organic-rich shale presents challenge in development of physics based predictive model for poromechanical properties. Owing to the direct economic impacts of shale rocks, there is an emerging necessity for developing multiscale models of poroelastic and damage behavior that bridge the coarse scale description of shale material to its structure and composition at subscale levels.
In the context of modeling poroelastic properties, multiscale model presented in (Ulm et al., 2005; Ulm and Abousleiman, 2006) offers a micromechanics based approach in which the homogenized properties at a coarser scale are result of applying effective medium approximations to a material volume element with fine scale representation. Utilizing this multiscale framework, the effect of texture, mechanical, thermal maturity and chemical composition on poroelastic behavior of shale are studied in some recent works (Abedi et al., 2016; Monfared and Ulm, 2016). In order to account for uncertainty emerging from inherent heterogeneity, Mashhadian et al. (2018) proposed a probabilistic model approach consisting of experimental characterization, micro-poromechanical modelling and uncertainty quantification and propagation. This development present promising framework to improve the predictive capability in modeling multiscale poroelastic properties of undamaged shale at macroscale. Quasi-brittle rock, such as shale, are susceptible to be weakened by various defects (e.g. microcracks, pores). The role of microcracks in the nonlinear mechanical behaviors of quasi-brittle rocks has been identified by many works (Simmons and Richter, 1976; Paterson, 1978). Based on standard Eshelby-based homogenization procedure Zhu et al. (2008) aims to describe the mechanical behavior of cracked solid by incorporating standard framework of irreversible thermodynamics. From then on various extension to the work has been done in several aspects (Zhu et al., 2011; Xie et al., 2011; Chen et al., 2014). However, most models developed so far have been essentially applied to dry materials. Few studies are conducted on micromechanical modeling of damage in saturated geo-materials. Dormieux and his coworkers were among founding researchers to have developed micromechanics of damage in saturated brittle media (Dormieux et al., 2006; Dormieux and Kondo, 2007). Along this line Xie et al. (2012), presents a poro-micromechanics based damage model which incorporates thermodynamics based damage criterion, and its application for saturated porous materials.
Giant reservoirs such as Lula (Santos Oil Basin, Brazil) and Ghawar (Saudi Arabia) have high-permeability intervals, known as super-k zones, associated with thin layers. Modeling these small-scale flow features in large-scale simulation models is difficult. Current methods are limited by high computational costs or simplifications that mismatch the representation of these features in simulation gridblocks. This work has two purposes: present an upscaling work flow to integrate highly laminated or interbedded reservoirs with thin, highly permeable layers in reservoir simulations through a combination of an explicit modeling of super-k layers using the Parsons (1966) formula and dual-medium flow models, and compare this method with two conventional upscaling approaches that are available in commercial software.
We use the benchmark model UNISIM-II-R (Correia et al. 2015a), a fine single-porosity grid dependent on field information from the Brazilian presalt and Ghawar oil fields, as the reference solution to compare the upscaling matching between the three methods. We compare oil recovery factor (ORF), water cut (WC), average reservoir pressure (RP), water front, and the time consumption for simulation. Our proposed Parson’s dual-medium (PDP) methodology achieved better upscaling matches with the reference model and had minimal time consumption compared with the representation of super-k layers through an implicit matrix modeling by single-porosity flow models (IMP) and through the explicit representation of super-k zones in the fracture system of dual-medium flow models (DFNDP).
Giant reservoirs such as Lula (Santos Oil Basin, Brazil) and Ghawar (Saudi Arabia) have high permeability intervals, known as super-k zones, associated with thin layers. Modeling these small-scale flow features in large-scale simulation models is complex. Current methods are limited by high computational costs or simplifications that mismatch the representation of these features in simulation grid blocks. This work has two purposes: (1) present an upscaling workflow to integrate highly laminated or inter-bedded reservoirs with thin, highly permeable layers in reservoir simulations through a combination of (a) an explicit modeling of super-k layers using Parsons (1966) formula and (b) dualmedium flow models, and (2) compare this method with two conventional upscaling approaches, available in commercial software. We use the benchmark model UNISIM-II-R, a fine single-porosity grid based on field information from the Brazilian Pre-salt and Ghawar oil fields, as the reference solution to compare the upscaling matching between the three methods. We compare; oil recovery, water cut, average reservoir pressure, waterfront, and the time consumption for simulation. Our proposed parsons dual-medium (PDP) methodology achieved better upscaling matches with the reference model and had minimal time consumption when compared with the representation of super-k layers through an implicit matrix modelling by single porosity flow models (IMP) and through the explicit representation of super-k zones in the fracture system of dual-medium flow models (DFNDP).
Steady-state two-phase relative permeability upscaling in synthetic and X-ray computerized tomography (CT) coal cores is performed with a three-dimensional (3D) reservoir simulator using an automated control procedure to drive a series of steady-state fractional flows. A clear understanding of relative permeability in coal is important for coalbed methane reservoir management from pore scale to sales point, as it is valuable for helping forecast production. Automation control enables greater continuity between physical corefloods and the numerical upscaling of the same coreflood procedure. Absolute permeability is computed for primitive synthetic core types using the reservoir simulator and is compared to an analytical formulation to validate the use of the simulator solution for core scale property determination. Relative permeability was computed for synthetic cores considering several scenarios: fracture geometry/ abundance (parallel vs. intersecting), rock-type matrix distribution (homogeneous vs. heterogeneous), ratios of matrix-to-fracture permeability (high vs. low), and injection rate conditions [capillary limit (CL) vs. viscous limit (VL)]. Additionally, injection rate conditions were evaluated in the upscaled relative permeability of an X-ray CT segmented composite coal core. Analysis of the upscaled relative permeability curves in the composite and synthetic cores illustrated the impact of each scenario on upscaling relative permeability and suggests that selected characteristics of unconventional cores can potentially be used to delineate parameter dependence in a manner similar to rock type volume fraction and ordering in conventional cores. The consistency of the developed upscaled results with previous studies confirms the applicability of automated process control in core scale multiphase upscaling using a commercial reservoir simulator at varied injection rates and upscaling conditions.
Wang, Lu (Lawrence Livermore National Laboratory) | Osei-Kuffuor, Daniel (Lawrence Livermore National Laboratory) | Falgout, Rob (Lawrence Livermore National Laboratory) | Mishev, Ilya (ExxonMobil) | Li, Jizhou (ExxonMobil)
Applications in geosciences, such as reservoir modeling, continue to grow in both size and complexity. Simulations are increasingly complex as they couple more physical phenomena over larger physical domains. As a result, the linear system that arises from the numerical solution of these problems can be challenging to solve by iterative methods. Simulators require advanced algebraic solvers that are robust enough to handle the anisotropies, heterogeneities and coupling between the physical variables, and scalable enough to handle large-scale solution on high performance parallel systems. Multigrid solvers are a class of iterative solvers that are scalable and efficient for solving linear systems that arise from large-scale applications. However, applications with multiple physical unknowns pose a challenge for standard multigrid techniques, particularly when the coupling between the unknowns is strong. In this paper, we present our efforts to develop a multigrid-preconditioned Krylov solver, where the preconditioner is based on the multigrid reduction framework. This preconditioner is designed to represent the coupling between the physical variables of the reservoir modeling equations and account for the underlying physics of the system. Two-stage preconditioners, such as the well-known constrained pressure residual (CPR) approach and its variants like CPR-AMG, have been commonly used in reservoir simulation applications. We discuss how these current solver strategies may be interpreted within the multigrid reduction framework to better understand the different variations. Finally, we present techniques for improving the MGR approach and present results on solver performance on examples from reservoir modeling.
A dynamic multilevel compositional solver (C-ADM) is introduced for fully(and sequentially-) implicit systems arising from compositional displacements in natural porous media. The fully (or sequential) implicit system is first described at a fine-scale resolution, where phases are allowed to consist of different components (based on thermodynamics equilibrium). In addition, heterogeneous capillary functions (defined based on Leverett's J-function) and gravitational effects are both considered, adding significantly to the nonlinear complexity of the processes. Given this complex fine-scale system for a heterogeneous reservoir, C-ADM defines a dynamic multilevel system, based on an error criterion, where the grid resolution is defined based on the physics of the process as well as geological complexities and location of wells. Once this multilevel grid is defined, sequences of prolongation and restriction operators are employed to obtain an accurate and efficient multilevel system. CADM allows for a general set of prolongation operators, e.g., constant, bilinear (or polynomial), and multiscale basis functions. The restriction operators, however, are constructed based on a mass-conservative finite-volume formulation at all levels. For several challenging test cases it is shown that C-ADM employs only a small fraction of the fine-scale grids to provide an accurate description of the process. C-ADM casts a promising approach in the application of dynamic grid refinement methods for real-field applications.