Interest in understanding the underlying physical mysteries that result in production from liquid rich shale (LRS) reservoirs is growing globally. Compositional modeling of fractured multiphase shale reservoirs is a complex process that is yet to be fully understood. Our quest in this work is to explore some of the key fundamentals of flow in nano-porous media. This study analyzes the impact of capillary pressure and critical property shifts due to the pore proximity effect on hydrocarbon production in multiple realistic scenarios of flow in shale reservoirs. While many past researches have attempted to show the significance of nano-forces in shale reservoirs, the actual impact that these have on large-scale hydrocarbon production is not fully understood. Response of wells in shale oil reservoirs is tied to both reservoir conditions and fluid system of interest. This work provides an assessment of the magnitude and the ultimate significance of the forces that arise in tight pores, using a fully compositional reservoir simulator. We studied the impact of capillary pressure amidst variations in pore size distributions, fluid composition, reservoir pressure conditions and the presence of fractures that are typical to shale reservoirs. In conventional compositional modeling phase behavior is assumed to be independent of capillary pressure. However, in the presence of nano-pores the magnitude of the differences in phase pressures can be very large. Here we present a systematic study of the influence of this phenomenon on recovery from shale reservoirs.
Our findings indicate that the impact of these nano-forces on hydrocarbon production is influenced by variations in shale reservoir/fluid properties. Shale reservoirs entail a great deal of geological uncertainty, and can contain a wide spectrum of hydrocarbon fluids that may range from low GOR black oil to dry gas. Further complications arise from uncertain pore size distribution and improperly characterized fracture networks.
Reservoir simulation models are used as an everyday tool for reservoir management. The number of grid blocks in a simulation model is generally much smaller than the number of the grid blocks in the geological model. The geological model is therefore regularly upscaled for building reasonably sized simulation models. Any upscaling causes a loss in detail and introduces errors. Understanding the nature of the errors that occur due to the upscaling step is important because it can provide a handle on the kind of upscaling to be performed and the optimum level of upscaling. Furthermore, understanding interaction of gridding and upscaling errors is essential for building reliable simulation models. In this paper we analyze errors introduced due to the upscaling step.
Firstly, the generality of the purely local single phase upscaling is considered and the errors introduced due to its use in complex multiphase flows are investigated. Three kinds of errors, namely total upscaling errors, discretization errors, and errors due to the loss of heterogeneity are defined and their behavior as a function of the level of upscaling is studied for different types of permeability distributions. The reasons for apparently low total upscaling errors at high levels of upscaling are investigated. The efficacy of the geostatistical tools (variograms, QQ plots) for understanding the geological structure of the upscaled models as compared to the reference model is tested for different cases.
Secondly, the uncertainty introduced in the flow results due to the upscaling errors is studied in conjunction with the uncertainty incorporated through the introduction of geological variability in the ensemble of geological models. The behavior of upscaling errors and geological variability as a function of level of upscaling is studied, and its possible impact on important reservoir management decisions, is analyzed. It is shown that the uncertainty models of cumulative oil profiles, obtained from flow results of highly upscaled models can contain significant uncertainties. These uncertainties are a result of upscaling errors and therefore cannot be considered to represent only geological uncertainty, which is captured by the introduction of geological variability in the ensemble of the geological models.
Geological models obtained through detailed reservoir characterization honor the structure and the petrophysics of the reservoir, but are rarely conditioned to the production data. Geological models are generally fine scale models preferably on the same scale as that of data acquisition, with the reservoir properties populated geostatistically. Reservoir simulation models make use of these geological models but the geological complexities of the fine scale models are often approximated to a certain extent by the upscaling of the fine scale model to a coarse scale model. These coarse scale models are coupled with the fluid and rock properties of the reservoir and other engineering data to obtain reservoir simulation models which are calibrated to field performance through history matching and then used to predict the performance of the reservoir under different operating conditions.
The level of detail incorporated in the reservoir simulation model depends not only on the time constraints and the ability of the available computer resources to efficiently handle the simulation model, but also on the study objectives. Number of grid blocks in a simulation model is an important parameter where the level of detail often available far exceeds the level of detail justified by the resource constraints and the study objectives. As a result of this, upscaling of the rock properties, mainly permeability, porosity and relative permeability, to reduce the number of grid blocks is a key step in a reservoir simulation study.
The paper deals with the upscaling of geological model wherein errors introduced due to the use of purely local single phase upscaling in complex multiphase flows is investigated. The error behavior as a function of the level of upscaling is investigated for isotropic and anisotropic permeability fields. The behavior of geological variability as a function of level of upscaling is studied in the light of its impact on the uncertainty.
Practical production optimization problems typically involve large, highly complex reservoir models, thousands of unknowns and many nonlinear constraints, which makes the numerical calculation of gradients for the optimization process impractical. This work explores a new algorithm for production optimization using optimal control theory. The approach is to use the underlying simulator as the forward model and its adjoint for the calculation of gradients. Direct coding of the adjoint model is, however, complex and time consuming, and the code is dependent on the forward model in the sense that it must be updated whenever the forward model is modified.
We investigate an adjoint procedure that avoids these limitations. For a fully implicit forward model and specific forms of the cost function and nonlinear constraints, all information necessary for the adjoint run is calculated and stored during the forward run itself. The adjoint run then requires only the appropriate assembling of this information to calculate the gradients. This makes the adjoint code essentially independent of the forward model and also leads to enhanced efficiency, as no calculations are repeated. Further, we present an efficient approach for handling nonlinear constraints that also allows us to readily apply commercial constrained optimization packages. The forward model used in this work is the General Purpose Research Simulator (GPRS), a highly flexible compositional/black oil research simulator developed at Stanford University.
Through two examples, we demonstrate that the linkage proposed here provides a practical strategy for optimal control within a general-purpose reservoir simulator. These examples illustrate production optimization with conventional wells as well as with smart wells, in which well segments can be controlled individually.
Most of the existing major oilfields are already at a mature stage, and the number of new significant discoveries per year is decreasing . In order to satisfy increasing worldwide demand for oil and gas, it is becoming increasingly important to produce existing fields as efficiently as possible, while simultaneously decreasing development and operating costs. Optimal control theory is one possible approach that can be deployed to address these difficult issues. The main benefit of the use of optimal control theory is its efficiency, which makes it suitable for application to real reservoirs simulated using large models, in contrast to many existing techniques.
The above-mentioned problem is essentially an optimization problem, wherein the objective is to maximize or minimize some cost function J(u) such as net present value (NPV) of the reservoir, sweep efficiency, cumulative oil production, etc. Here, u is a set of controls including, for example, well rates and bottom hole pressures (BHP), which can be manipulated in order to achieve the optimum. In other words, u is anything that can be controlled. It should be understood that the optimization process results in control of future performance to maximize or minimize J(u), and thus the process of optimization cannot be performed on the real reservoir, but must be carried out on some approximate model. This approximate model is usually the simulation model of the reservoir. This simulation model is a dynamic system that relates the controls u to the cost function J(u).
Consider for example the simple schematic of a reservoir shown in Fig. 1, where the cost function is cumulative oil production and the control is the injection rate. Changing the injection rate changes the dynamic states of the system (pressures, saturations), which changes the oil production rate of the producer, which in turn impacts the cost function. Thus, the controls u are related to J(u) through the dynamic system. The dynamic system can also be thought of as a set of constraints that determine the dynamic state of the system given a set of controls. Further, the controls u themselves may be subject to other constraints that dictate the feasible or admissible values of the controls, such as surface facility constraints or fracture pressure limits. It is these additional constraints that in many cases complicate the problem and the solution process.
Drift-flux models represent multiphase flow in wellbores or pipes in terms of a number of empirically determined parameters. Because of the lack of data for two and three-phase flow in large-diameter inclined pipes, existing parameters are commonly based on small-diameter pipe experiments, which can lead to significant errors when the models are applied to wellbore flows. In this work, we use recent large-diameter experimental data for the determination of drift-flux parameters for oil-water-gas flow. The parameters are computed through application of an optimization procedure. It is shown that gas holdup in three-phase systems can be estimated using a two-phase flow model by viewing the system as an effective gas-liquid system, with oil and water comprising the ‘liquid' phase. This approach is, however, generally inaccurate for the determination of oil and water holdups, in which case the effect of gas must be taken into account. Specifically, for pipe inclinations away from horizontal, even small amounts of gas can act to eliminate the slip between oil and water. As the pipe deviation approaches horizontal, however, oil-water slip persists, even in the presence of gas. We develop and apply a unified two and three-phase flow model to capture this gas effect. The new model is shown to provide much more accurate predictions for oil and water holdups in three-phase systems than were achievable with previous models.
A simplified discrete-fracture model suitable for use with general-purpose reservoir simulators is presented. The model handles both 2- and 3D systems and includes fracture-fracture, matrix-fracture, and matrix-matrix connections. The formulation applies an unstructured control volume finite-difference technique with a two-point flux approximation. The implementation is generally compatible with any simulator that represents grid connections by a connectivity list. A specialized treatment based on a "star-delta" transformation is introduced to eliminate control volumes at fracture intersections. These control volumes would otherwise act to reduce numerical stability and timestep size. The performance of the method is demonstrated for several oil/water flow cases including a simple 2D system, a more complex 3D fracture network, a localized fractured region with strong capillary pressure effects, and a model of a strike-slip fault zone. The discrete-fracture model is shown to provide results in close agreement with those of a reference finite-difference simulator in cases in which direct comparisons are possible.
Shi, H. (Stanford University) | Holmes, J.A. (Schlumberger GeoQuest) | Durlofsky, L.J. (ChevronTexaco EPTC) | Aziz, K. (Stanford University) | Diaz, L.R. (Stanford University) | Alkaya, B. (Stanford University) | Oddie, G. (Schlumberger Cambridge Research)
Drift-flux modeling techniques are commonly used to represent two and three-phase flow in pipes and wellbores. Unlike mechanistic models, drift-flux models are continuous, differentiable and relatively fast to compute, so they are well suited for use in wellbore flow models within reservoir simulators. Drift-flux models require a number of empirical parameters. Most of the parameters used in current simulators were determined from experiments in small diameter (2 inch or less) pipes. These parameters may not be directly applicable to flow in wellbores or surface facilities, however, as the flow mechanisms in small pipes can differ qualitatively from those in large pipes. In order to evaluate and extend current drift-flux models, an extensive experimental program was initiated. The experiments entailed measurement of water-gas, oil-water and oil-water-gas flows in a 15 cm diameter, 11 m long plexiglass pipe at 8 deviations ranging from vertical to slightly downward. In this paper, these experimental data are used to determine drift-flux parameters for steady state two-phase flows of water-gas and oil-water in large-diameter pipes at inclinations ranging from vertical to near-horizontal. The parameters are determined using an optimization technique that minimizes the difference between experimental and model predictions for holdup. It is shown that the optimized parameters provide considerably better agreement with the experimental data than do the existing default parameters.
A simplified discrete fracture model suitable for use with general purpose reservoir simulators is presented. The model handles both two and three-dimensional systems and includes fracture-fracture, matrix-fracture and matrix-matrix connections. The formulation applies an unstructured control volume finite difference technique with a two-point flux approximation. The implementation is generally compatible with any simulator that represents grid connections via a connectivity list. A specialized treatment based on a "star-delta" transformation is introduced to eliminate control volumes at fracture intersections. These control volumes would otherwise act to reduce numerical stability and time step size. The performance of the method is demonstrated for several example cases including a simple twodimensional system, a more complex three-dimensional fracture network, and a model of a strike-slip fault zone. The discrete fracture model is shown to provide results in close agreement with those of a reference finite difference simulator in cases where direct comparisons are possible.
Flow through fractured porous media is typically simulated using dual-porosity models. This approach, although very effi- cient, suffers from some important limitations. One drawback is that the method cannot be applied to disconnected fractured media. In addition, dual-porosity models are not well suited for the modeling of a small number of large-scale fractures, which may dominate the flow. Another shortcoming is the difficulty in accurately evaluating the transfer function between the matrix and the fractures. For these reasons, discrete fracture models, in which the fractures are represented individually, are beginning to find application in reservoir simulation. These models can be used both as stand-alone tools as well as for the evaluation of transfer functions for dual-porosity models. Such models can also be used in combination with the dual-porosity approach.
To accurately capture the complexity of a fractured porous media it is usually necessary to use an unstructured discretization scheme. There are, however, some effective procedures based on structured discretization approaches. For example, Lee et al.1 presented a hierarchical modeling of flow in fractured formations. In this approach the small fractures were represented by their effective properties and the large-scale fractures were modeled explicitly. In the case of unstructured discretizations there are two main approaches - finite element and finite volume (or control volume finite difference) methods. Baca et al.2 were among the early authors to propose a two-dimensional finite element model for single phase flow with heat and solute transport. In a more recent paper Juanes et al.3 presented a general finite element formulation for two- and three-dimensional single-phase flow in fractured porous media.
There has been some work on the extension of the finite element method to handle multiphase flow. For example, Kim and Deo4 and Karimi-Fard and Firoozabadi5 presented extensions of the work of Baca et al.2 for two-phase flow. They modeled the fractures and the matrix in a two-dimensional configuration with the effects of capillary pressure included. The two media (matrix and fractures) were coupled using a superposition approach. This entails discretizing the matrix and fractures separately and then adding their contributions to obtain the overall flow equations.
Coupled geomechanical-fluid flow models are needed to account for rock deformation resulting from flow-induced pressure changes in stress sensitive reservoirs. There are, however, issues of numerical stability that must be addressed before these coupled models can be used reliably. Specifically, it is known that standard procedures can lead to pressure oscillations as a result of the violation of the Babuska-Brezzi (B-B) condition, which requires unequal order interpolation of the displacement and pore-pressure variables. In this paper, a number of different types of coupled models are considered. A novel finite element method is developed to circumvent the B-B condition. The method applies a stabilized finite element technique to solve the force balance and pressure equations along with a finite volume method to solve the remaining component mass balance equations. All of the equations are solved in a fully coupled fashion. This method is compared with fully coupled and iteratively coupled models developed using non-stabilized finite elements for the force balance with finite volume methods for all of the component mass balance equations. These comparisons demonstrate that all methods perform reliably on homogeneous reservoirs over long time scales. The stabilized method is shown to provide improved stability at early times and for reservoirs with very low permeability barriers.
Geomechanics can play an important role in oil recovery and oil field operations. In particular, geomechanics is crucial in such problems as production-induced compaction and subsidence, borehole stability, and hydraulic fracturing. In these areas, a key mechanism is the rock deformation due to changing in situ stress. Several investigations1,2,3 have shown that models that do not couple flow and geomechanics may give inaccurate predictions. For this reason, recent attention has been focused on modeling of the coupled system1,2,3,4,5.
Here we present an accurate and stable numerical technique for coupled problems. In order to build robust and accurate numerical models, it is important to first consider the accuracy and stability properties of existing methods. Vermeer and Verruijt6 presented a stability condition for a simple consolidation problem. This condition required that time steps be larger than a certain value, otherwise spatial oscillations would occur. Such a requirement may seem counterintuitive, as we generally expect to observe improved accuracy as the time step is decreased. Within the context of coupled soil mechanics - fluid flow problems, Zienkiewicz et al.7 showed that the problem could be understood in terms of the Babuska-Brezzi stability condition (B-B condition)8,9, which requires nonuniform (i.e., different order) interpolation of the displacement and pore-pressure variables to ensure stability. Murad and Loula10,11 showed that an incompressible Stokes flow system at early times led to spatial instability if the B-B condition was not satisfied.
Hughes et al.12 developed stabilized methods to circumvent the B-B restrictions for incompressible Stokes flow problems. These methods are within the category of weighted residual methods. Weighted residual methods modify the bilinear form associated with the problem so that improved numerical stability is achieved without compromising consistency. Zienkiewicz andWu13 showed similar stabilization by modifying the form of the time stepping in transient problems.
In this work, we adopt the general approach of stabilized methods and develop stabilized finite element schemes to control the numerical instability arising in consolidation problems without losing consistency. Using this approach, any combination for interpolations of pressure and displacement can be used. We demonstrate both the instability that can be observed using standard methods as well as the improved performance of our stabilized method.
A well model is required to relate the well rate to the well pressure and the wellblock pressure while modeling wells in reservoir simulators. The well index is calculated automatically by simulators for conventional wells, but it is generally calculated by the user and supplied in the input while modeling nonconventional wells when the default procedure available in commercial simulators is not adequate.
In this paper we describe a new semi-analytical solution for horizontal wells with multiple fractures. The fractures can be rotated at any horizontal angle to the well, and they need not fully penetrate the formation in the vertical direction. This solution for hydraulically fractured wells is obtained by applying Fourier analysis to a 2D solution; therefore, this solution is easy to obtain when the 2D solution is available. The analytical solution provides the well pressure that can be combined with a numerically computed gridblock pressure to obtain the well index (WI). Results of our rigorous solution are compared with empirical approaches currently in use for calculating the well index and the productivity index. Examples are given when horizontal wells with fractures are modeled using our new approach and conventional methods. These results demonstrate the need for computing the correct WI.
The technology of fracturing horizontal wells is being widely used by the petroleum industry. Therefore one will ask the questions: Are the currently available models adequate for predicting performance of horizontal wells with fractures? If not, how can we develop more appropriate models?
For the modeling of wells in reservoir simulators, a well model is required to relate the well flow rate to the well pressure and the wellblock pressure. Peaceman1,2 established a mathematical relationship between the wellblock pressure and wellbore pressure for a fully penetrating vertical well under the condition of 2D flow. This relationship is valid for isolated wells under steady-state or pseudosteady conditions. Peaceman's well model is usually the default in general reservoir simulators. Because of the variety of nonconventional wells used by industry, such as horizontal and fractured wells, and because of the variety of grid systems employed, much research is being conducted on improving simulator accuracy. Some of the more important works are summarized here.
Babu et al.3 extended Peaceman's work to the case of a uniform flux horizontal well in a slab-like drainage area. Penmatcha et al.4,5 have extended Babu and Odeh's model to the case of infinite well conductivity. Jasti et al.6 present results to increase simulator accuracy when solving the problem of a system of partially penetrating wells having arbitrary trajectories. Ding7 used the concept of layer potentials and applied transmissibility adjustments. Maizeret8 developed an analytical solution to investigate well indices for nonconventional wells. In 1998, Wan et al.9 compared several well models, including explicit well models, and showed the importance of correct well indices.
Our focus here is on horizontal wells with single or multiple fractures. Basically, there are three ways to model hydraulic fractures:
Refine the grids, so that the true geometry of the fractures is modeled. In this method, fractures are represented by very fine gridblocks; properties of the fractures, such as permeability, are assigned to those fracture blocks. Schulte10 evaluated the transient pressure behavior of a fractured vertical well in this way. However, the size of the fracture blocks was so small that numerical instabilities were encountered. To overcome this problem, Roberts et al.11 increased the fracture grid width and decreased the fracture permeability.
Modify the transmissibilities of the blocks that contain the fractures. Work has been done in independently generating fracture transmissibilities with the transmissibilities of reservoir blocks adjacent to the fracture plane.12,13 Therefore, the productivity of the fracture is included while the negligible volume of the fracture is excluded. These authors also apply various coupling methods to deal with the interfaces between the fracture model and the reservoir model. Hegre14 computed correct transmissibilities for a coarse grid, which gives the same long-term pressure behavior as a fine-grid model. The transmissibility adjustments are based on results from a model, which is fine gridded in the near wellbore/ fracture area with explicit modeling of the fractures.
Modify the effective wellbore radius or well index. Cinco- Ley and Samaniego15 showed that the effective wellbore radius is one-fourth of the total fracture length for an infinite conductivity vertical fracture, so this effective wellbore radius can be used instead of rw in the productivity index formulas. Hegre and Larsen16 showed that the effective wellbore radius for a horizontal well intercepted by either a longitudinal fracture or multiple transverse fractures can be analytically determined. They used the rew calculated in this manner in the well index formula (see Eq. 15) instead of rw.
The approach of modifying the effective wellbore radius is used in this study. A semi-analytical approach is used to calculate the well index.6,8 First we develop the analytical solution for a horizontal well with multiple fractures in a brick-shaped drainage volume. With the wellbore pressure from this analytical solution and the wellblock pressure from a single-phase numerical simulator, the effective wellbore radius is calculated.
Complex well configurations coupled with intelligent completions offer great potential for the efficient production of oil reservoirs. One application of this technology is the use of surface adjustable downhole chokes. This allows wells to be divided into separate inflow zones, with each zone produced under a different pressure drawdown, but with production commingled. Chokes can be set to provide a more uniform inflow profile, which acts to delay the breakthrough of water and gas. The practical benefits of this technology have been demonstrated in several North Sea installations.
Despite the large technical and economic risks associated with these complex wells, little work has been presented on the detailed modeling of their performance. While it is possible to apply existing finite difference reservoir simulators, such models can be time consuming to build and the accuracy of the results depends on the grid and the well model. Here, we apply semi-analytical solution methods (based on Green's functions), appropriate for modeling the performance of non-conventional wells operating under single phase flow conditions, to these complex well configurations. The approach entails a fully coupled nonlinear formulation that accounts for pressure drops in the annulus, tubing and downhole chokes.
Numerical results for a variety of cases are presented. The high degree of accuracy of the semi-analytical approach is first demonstrated through comparison to results from established methods for model problems. Then, using a geostatistical description of a highly heterogeneous fluvial reservoir, we demonstrate how the method can be used to model the performance of wells with intelligent completions. We show how the technique can be applied to determine downhole choke settings that provide optimal inflow profiles.
The use of horizontal and deviated wells has had a dramatic impact on the industry over the last decade. These wells greatly increase reservoir exposure and can increase the production from a single well significantly. However, horizontal well performance can be impacted significantly by reservoir heterogeneity. Fluvial reservoirs, for example, typically display very large differences in permeability depending on whether the rock type is a channel sand or not, with the result that most of the inflow into a horizontal well might originate from a small proportion of the perforated length. In layered or faulted reservoirs additional complications can arise; e.g., the formation permeability and pressure may vary so much that a single well completion might not be feasible.
These complications can be addressed through the use of complex well configurations coupled with so-called "intelligent completions." A key aspect of this technology is the use of surface adjustable downhole chokes, which allow for the isolation of individual reservoirs or zones within a single reservoir. The downhole system makes use of packers that are set between the liner and tubing. It is then possible to control the inflow from different zones along the wellbore through use of downhole chokes, across which a pressure drop can be imposed. This type of completion has been successfully installed and operated in several field applications.1
Most intelligent completions are very expensive, and it is important to be able to model and predict the performance of these complex systems. Finite difference reservoir simulators can be used to provide an in-depth analysis of the problem,2 though there are several complications with this type of simulation. Specifically, finite difference models have large data requirements, especially when introducing complex well configurations. As a result, the tool might not be appropriate for quick screening studies. Further, because the reservoir is discretized, well connection factors (or well indexes), linking the well to the discrete model, are required. Default connection factors are typically based on the traditional Peaceman expression,3,4 developed under the assumption of two dimensional flow in a homogeneous reservoir. These well indexes may not be valid for complex wells with intelligent completions operating in heterogeneous reservoirs.