Search

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.

In production optimization, we seek to determine the well settings (bottomhole pressures, flow rates) that maximize an objective function such as net present value. In this paper we introduce and apply a new approximate dynamic programming (ADP) algorithm for this optimization problem. ADP aims to approximate the global optimum using limited computational resources via a systematic set of procedures that approximate exact dynamic programming algorithms. The method is able to satisfy general constraints such as maximum watercut and maximum liquid production rate in addition to bound constraints. ADP has been used in many application areas, but it does not appear to have been implemented previously for production optimization. The ADP algorithm is applied to two-dimensional problems involving primary production and water injection. We demonstrate that the algorithm is able to provide clear improvement in the objective function compared to baseline strategies. It is also observed that, in cases where the global optimum is known (or surmised), ADP provides a result within 1-2% of the global optimum. Thus the ADP procedure may be appropriate for practical production optimization problems.

Advanced wells are being used more and more frequently to improve the economics of oil field operations. Such wells can contact larger regions of the reservoir than is possible with standard wells. Because these wells can be very long, pressure drop along the well and temperature effects (in thermal problems) can became quite important. Modeling of advanced wells where thermal effects need to be included is a difficult task, and available models are often limited in their capabilities or are inefficient. A suitable model should be able to handle complicated well trajectories, heat transfer to and from the well, the hold-up of fluids due to slip between phases, and pressure change due to friction. This paper presents a comprehensive multisegmented well model that includes these features. The model has been implemented in Stanford's General Purpose Research Simulator (GPRS). It calculates both pressure and temperature profiles along the well for standard wells and wells with complicated trajectories. One of the advantages of our model is that both homogeneous (no slip) and drift-flux flow models can be used to determine phase distribution within the well. The hydrocarbon fluids can be described either through a black oil representation or using a general compositional model with an arbitrary number of components. The proposed multisegmented well model was tested with a number of challenging and realistic examples. Our results show that the model is stable and in most cases converges in only a few iterations even with reasonably large time steps. The capabilities of the model are demonstrated using three examples. The first example involves a vertical production well with the fluid represented in terms of six hydrocarbon components and water. For this case we demonstrate the importance of slip between phases. The second example is for a thermal multilateral well with two branches and a fluid with three hydrocarbon components. The third example shows our ability to model phase appearance and disappearance in the wellbore and demonstrates the impact of temperature on oil and gas rates during production.

As many fields around the world are reaching maturity, the need to develop new tools that allow reservoir engineers to optimize reservoir performance is becoming more urgent. One of the more challenging and important problems along these lines is the well placement optimization problem. In this problem, there are many variables to consider: geological variables like reservoir architecture, permeability and porosity distributions, and fluid contacts; production variables, such as well placement, well number, well type, and production rate; and economic variables like fluid prices and drilling costs. Furthermore, availability of complex well types, such as multilateral wells (MLWs) and maximum reservoir contact (MRC) wells, aggravate this challenge. All these variables, together with reservoir geological uncertainty, make the determination of an optimum development plan for a given field difficult. The objective of this work was to employ an optimization technique that can efficiently address the aforementioned challenges. Based on the success and versatility of Genetic Algorithms (GAs) in problems of high complexity with high dimensionality and nonlinearity, it is used here as the main optimization engine. Both binary GA (bGA) and continuous GA (cGA) were tested in the optimization of well location and design in terms of well type, number of laterals, and well and lateral trajectories in a channelized synthetic model. Both GA variants showed significant improvement over initial solutions but comparisons between the two types showed that the cGA was more robust for the problem under consideration. The cGA was, thereafter, applied to a real field located in the Middle East to investigate its robustness in optimizing well location and design in more complex reservoir models. The model is an upscaled version for an offshore carbonate reservoir, which is mildly heterogeneous with low and high permeability areas scattered over the field. After choosing the optimization technique to achieve our objective, considerable work was performed to study the sensitivity of the different algorithm parameters on converged solutions. Then, multiple optimization runs were performed to obtain a sound development plan for this field. An attempt was made to quantify how solutions were affected by some of the assumptions and preconditioning steps taken during optimization. Finally, an optimization ran was performed on the fine model using optimized solutions from the coarse model. Results showed that the optimum well configuration for the reservoir model at hand can contain five or more laterals; which shows potential for drilling MRC wells. Other studies comparing results from the fine and coarse reservoir models revealed that the best solutions are different between the two models. In general, solutions from different runs had different well designs due to the stochastic nature of the algorithm but some guidance about preferred well locations could be obtained through this process

Thermal recovery processes are widely used for heavy oil production and are under investigation for other unconventional resources. The management and optimization of these processes require accurate representations of relevant physical phenomena in simulation models. An important aspect of the simulation capability is the modeling of flow in the wellbore and the coupling of wellbore and reservoir flows. In previous work, we developed and applied a black-oil thermal multiphase wellbore model and linked it to Stanford's General Purpose Research Simulator (GPRS). In this paper, we extend this formulation to treat thermal compositional systems. GPRS already contains a thermal compositional capability for reservoir flow; here we develop and test a thermal compositional wellbore model. This wellbore model includes an energy conservation equation, mass conservation equations for each component, and a general pressure drop relationship. The multiphase wellbore flow is represented using a drift-flux model, which includes slip between the three phases. The model determines the temperature, pressure, mixture flow velocity and component fractions as functions of time and axial position along the well. We apply the coupled thermal compositional wellbore-reservoir model to the simulation of several cases involving thermal effects. These include a verification example, in which we compare results from the compositional formulation to those of a black-oil model, and an example with seven components. Different well geometries (e.g., dual-lateral and deviated with varying inclination angle) are considered in these applications.

Field development optimization is a computationally intensive task due to the large number of reservoir simulation runs required. These simulations can be expensive, especially for large and complex reservoir models. Proxies can be used to efficiently estimate the objective function value for new scenarios and can act to reduce the number of simulations required. Thus they can be very useful for speeding up field development optimization.In this paper a procedure that combines an optimization algorithm (in this case a genetic algorithm or GA) and a new statistical proxy is described. The statistical proxy has the following key elements. First, a new selection procedure called individual-based selection is applied to decide which individuals (scenarios) are to be simulated. Second, the new approach uses multiple proxies for optimization problems involving multiple reservoir models, which are needed to account for geological uncertainty. Third, the statistical proxy is modified to work efficiently in distributed computing environments. Finally, the proxy procedure is successfully incorporated into an existing general field development optimization package (Williams et al., 2004; Litvak et al., 2007a).In the individual-based selection method, for each scenario the proxy estimate of the objective function is compared to a threshold. If the estimate exceeds the threshold, then the case is simulated (otherwise it is not simulated). The threshold corresponds to a specified percentile of the cumulative distribution function constructed from previously simulated cases and therefore changes during the course of the optimization. In cases with multiple reservoir models, each model has its own corresponding proxy. This eliminates the problem of duplicate objective function estimates for different reservoir models, which may occur with previous proxy-based methods. The individual-based selection method is shown to perform better for a particular example than the population-based method published previously.The overall procedure is applied to the optimization of infill drilling where we maximize the incremental net present value (NPV) by optimizing new well locations, well type and rig schedule, subject to field development constraints. We demonstrate the capabilities of the proxy using synthetic reservoir models and a real field in the Gulf of Mexico. In the first example, two optimization cases are considered, corresponding to the use of single and multiple reservoir models. In the case with one reservoir model, the hybrid procedure found the same field development scenario compared to GA only, and required 85% fewer simulations. In the case with multiple reservoir models, the hybrid procedure found a slightly different field development scenario than the pure GA approach, though the NPV from the hybrid procedure was within 1% of that using only GA. The hybrid approach, however, required 91% fewer simulations for this case. In the field application, a better field development scenario with 45% fewer simulations was found using the hybrid algorithm (GA and proxy) compared to using only GA. These examples clearly demonstrate the effectiveness of the statistical proxy procedure for accelerating field development optimization.

The detailed interactions between the reservoir and the wellbore are especially important in thermal processes such as steam flooding and in situ upgrading. These linkages must therefore be captured in thermal simulations. Although fully-coupled thermal wellbore-reservoir flow simulators have been developed, the implementation of the thermal well model is somewhat complicated and the simulations are computationally demanding. In this paper, we present a semianalytical treatment that enables a fairly straightforward extension of existing isothermal wellbore flow models to the nonisothermal case. The procedure entails the use of analytical solutions for wellbore temperature applied in conjunction with numerical solutions of the reservoir mass and energy balance equations coupled with wellbore mass and momentum balance equations. The approach thus enables a degree of decoupling between the wellbore flow and energy problems. We proceed by first presenting analytical solutions for wellbore temperature, developed under various assumptions (some of these solutions have been obtained previously). We then describe the use of one of these solutions, which allows for general variation of in situ phase fraction and other properties along the wellbore, within the semianalytical context. The implementation of the overall method into a general purpose research simulator is also described. Results are presented for several cases involving multiphase flow in monobore and multilateral wells. Close agreement with reference solutions, obtained from a fully-coupled thermal wellbore-reservoir model, is demonstrated for all of the examples. The semianalytical treatment is additionally shown to provide comparable or improved computational efficiency relative to the fully-coupled model. The overall procedure is therefore very well suited for use in general thermal reservoir simulation.

Efficient history matching of geologically complex reservoirs is important in many applications, but it is central in closed-loop reservoir modeling, in which real-time model updating is required. Within the context of closed-loop reservoir modeling, the two approaches receiving the most attention to date are ensemble Kalman filtering and gradient-based methods using Karhunen-Loeve representations (eigen-decomposition) of the permeability field. Both of these procedures are technically appropriate only for random fields (e.g., permeability) characterized by two-point geostatistics (multi-Gaussian random fields). Realistic systems are much better described by multipoint geostatistics, which is capable of representing key geological structures such as channels. History matching algorithms that are able to reproduce realistic geology provide enhanced predictive capacity and are therefore more suitable for use with field optimization. In this work, we apply a new parameterization, referred to as a kernel principal component analysis (kernel PCA or KPCA) representation, to model permeability fields characterized by multipoint geostatistics. Kernel PCA enables preserving arbitrarily high order statistics of random fields, thereby providing the capability to reproduce complex geology. The KPCA representation is then combined with an efficient gradient-based history matching technique. The linkage of KPCA for modeling geology with gradient-based history matching is very natural as the KPCA representation is differentiable and gradients with respect to geological parameters can be readily computed. The overall procedure is then applied to several example cases, including synthetic models and a model of a real reservoir. The approach is shown to better reproduce complex geology, which leads to improved history matches and better predictions, while retaining reasonable computational requirements.

Abstract We develop and analyze an algorithmic framework and library of parallel automatically differentiable data-types written in the C++ programming language. In addition to representing variables and discretized expressions on a simulation grid, the data-types hide algorithms employed behind the scenes capable of automatically computing the sparse analytical Jacobian. Using the library, reservoir simulators can be developed rapidly by simply writing the residual equations, and without any hand differentiation, loops, or any other low-level constructs. The key challenge addressed is a technical development to enable this to run fast. Faster than if any of several existing automatic differentiation packages were used, faster than any purely Object Oriented implementation, and at least as fast as if a development team implemented hand-optimized residuals, analytical derivatives, and Jacobian assembly routines. This is achieved by extending a generic programming technique specific to the C++ programming language known as Expression Templates. The extension enables its application to handle nested abstract types including sparse vector composites. Consequently, the library implementation can handle the complexity of any reservoir simulation model with optimal use of microprocessor hardware. We demonstrate the ease of use of the library by presenting an efficient and rapidly developed fully-implicit simulator. Performance analysis is performed with varying grid size and is bench-marked with structured implementations. The tests are run with widely used compilers on standard microprocessors. Introduction Next generation simulators aim to simultaneously (1) leverage current technology to solve today's problems most efficiently, and (2) be easily extendable to leverage tomorrow's technology to solve tomorrow's problems most efficiently. Towards this, the General Purpose paradigm makes conceptual abstraction of equations and unknowns to achieve extendability to a wide range of formulations and recovery processes 1,2,3,4. Another kind of conceptual abstraction is on numerical formulations and solution processes, and aims to address extendability to evolving computer technologies and solution algorithms. Motivation. Next generation simulator designs translate both conceptual abstractions into appropriate software-engineering abstractions. In this fashion, typical designs are made up of modules in a multi-layered framework 5,6,7,8. Each layer consists of modules with well defined interfaces that implement a certain level of functionality. Higher-layer modules implement increasingly complex functionality by using lower-layer modules. Apparent advantages of taller hierarchies are greater flexibility, extendability, and reliability.

This work discusses a closed-loop approach for efficient real-time production optimization that consists of three key elements – adjoint models for efficient parameter and control gradient calculation, polynomial chaos expansions for efficient uncertainty propagation, and Karhunen-Loeve (K-L) expansions and Bayesian inversion theory for efficient real-time model updating (history-matching). The control gradients provided by the adjoint solution are used by a gradient-based optimization algorithm to determine optimal control settings, while the parameter gradients are used for model updating. Polynomial chaos expansions provide optimal encapsulation of information contained in the input random fields and output random variables. This approach allows the forward model to be used as a black box but is much faster than standard Monte Carlo techniques. The K-L representation allows for the direct application of adjoint techniques for history matching while assuring that the two-point geostatistics of the reservoir description are maintained. The benefits and efficiency of the overall closed-loop approach are demonstrated through real-time optimizations of net present value (NPV) for synthetic reservoirs under waterflood subject to production constraints and uncertain reservoir description. The closed-loop procedure is shown to provide a substantial improvement in NPV over the base case, and the results are seen to be very close to those obtained when the reservoir description is known apriori.

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.

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.

Summary 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 in-situ gas volume fraction 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 constituting 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.

This paper discusses new techniques for the modeling and simulation of naturally fractured reservoirs with dual porosity models. Most of the existing dual porosity models idealize matrix-fracture interaction by assuming orthogonal fracture systems (parallelepiped matrix blocks) and pseudo-steady state flow. More importantly, a direct generalization of single-phase flow equations is used to model multi-phase flow, which can lead to significant inaccuracies in multiphase flow behavior predictions. In this work, many of these existing limitations are removed in order to arrive at a transfer function more representative of real reservoirs. Firstly, combining the differential form of the single-phase transfer function with analytical solutions of the pressure diffusion equation, an analytical form for a shape factor for transient pressure diffusion is derived to corroborate its time dependence. Further, a pseudo-steady shape factor for rhombic fracture systems is also derived and its effect on matrix-fracture mass transfer demonstrated. Finally, a general numerical technique to calculate the shape factor for any arbitrary shape of the matrix block (i.e. non-orthogonal fractures) is proposed. This technique also accounts for both transient and pseudo-steady state pressure behavior. The results were verified against fine-grid single porosity models and were found to be in excellent agreement. Secondly, it is shown that the current form of the transfer function used in reservoir simulators does not fully account for the main mechanisms governing multiphase flow. A complete definition of the differential form of the transfer function for two-phase flow is derived and combined with the governing equations for pressure and saturation diffusion to arrive at a modified form of the transfer function for two-phase flow. The new transfer function accurately takes into account pressure diffusion (fluid expansion) and saturation diffusion (imbibition), which are the two main mechanisms driving multiphase matrix-fracture mass transfer. New shape factors for saturation diffusion are defined. It is shown that the prediction of wetting phase imbibition using the current form of the transfer function can be quite inaccurate, which might have significant consequences from the perspective of reservoir management. Fine grid single porosity models are used to verify the validity of the new transfer function. The results from single block dual porosity models and the corresponding single porosity fine grid models were in good agreement.

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.

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 Babuška-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.

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 two-dimensional 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.

Downhole inflow control devices allow for the flexible operation of nonconventional wells. They provide the ability to independently control each branch of a multilateral well and can therefore be used to maximize oil production or to minimize the production of unwanted water or gas. In this paper, a method for the optimization of the operation of smart wells; i.e., wells containing downhole sensors and inflow control devices, is presented and applied. The method entails the use of a conjugate gradient optimization technique applied in conjunction with a commercial reservoir simulator that contains a detailed multisegment well model. The overall optimization technique is applied to several example problems involving different types of wells and geological models as well as multiple geostatistical realizations. The improvement in predicted performance using inflow control devices, which is as high as 65% in one case, is demonstrated for all of the examples considered. There is, however, significant variation in the level of improvement attainable using these devices, so sophisticated decision-making techniques may be required when considering their use in practice.

The current trend in industry is to use finer and finer grid with larger and larger number of components. This requires the development of faster simulation techniques than those currently available in commercial simulators. The fully implicit (FIM) model, the implicit pressure and explicit saturations (IMPES) model and the adaptive implicit (AIM) model are the traditional approaches. The IMPSAT (implicit pressure and saturations and explicit component mole fractions) model has been proposed in the literatures, but so far it is not considered a viable approach. This may be a result of lack of suitable stability analysis and improper implementation. Our analysis shows that the IMPSAT model is significantly more stable than the IMPES model, and in many cases substantially less expensive than the FIM model. The IMPSAT model becomes particularly attractive as the number of components becomes large. IMPSAT can also be used in an adaptive implicit approach to form various IMPSAT based adaptive implicit (IMPSAT-AIM) models. A new General Purpose Research Simulator (GPRS) using a General Formulation Approach has been developed and tested. The General Formulation Approach used in GPRS can be used to obtain any type of model. We can easily form a model by selecting any reasonable set of primary variables and any level of implicitness from IMPES to FIM. One important subset of models that have been developed from this new approach is the IMPSAT model and the IMPSAT based AIM models. The performance of different compositional models has been evaluated using a wide range of problems. We found that the IMPSAT model and the IMPSAT-AIM models have very good performance, and they are more efficient than the traditional models (FIM, IMPES and AIM). The IMPSAT model is generally over 50% faster than the IMPES model due to its improved stability. For problems that are not very hard, the IMPSAT model is cheaper to run than the FIM model since it reduces the number of unknowns that have to be solved implicitly. The IMPSAT-AIM models are more flexible, more stable and with appropriate linear solvers have the potential to be less expensive in terms of CPU time than the traditional AIM model.

For some petroleum fields, optimization of production operations can be a major factor for increasing production rates and reducing production costs. While for single wells or other small systems simple nodal analysis can be adequate, large complex systems demand a much more sophisticated approach. Many mature fields are produced by gas-lift under multiple constraints imposed by the field handling capacity of the system. In this paper we present an optimization technique for allocating production rates and lift-gas rates to wells of large fields subject to multiple flow rate and pressure constraints. The well rate and lift-gas rate allocation problem has been addressed in the literature1-13. However, existing methods are either inefficient or make significant simplifications. This often leads to suboptimal operations. This paper proposes a new formulation of the problem that is able to handle flow interactions among wells and can be applied to a variety of problems of varying complexities. We show that the proper formulation of the optimization problem is important in the practical use of modern optimization techniques. Once formulated, the optimization problem is solved by a sequential quadratic programming algorithm. Our results show that the procedure developed in this paper is capable of handling complex oil production problems.

Summary The determination of the optimal type, location, and trajectory of a nonconventional well is extremely challenging. The problem is more complicated than other well optimization problems because of the wide variety of possible well types (i.e., number, location, and orientation of laterals) that must be considered. In this paper, a general methodology for the optimization of nonconventional wells is presented. The optimization procedure entails a Genetic Algorithm (GA) applied in conjunction with several acceleration routines that include an artificial neural network, a hill climber, and a near-well upscaling technique. The overall methodology is then applied to a number of problems involving different reservoir types and fluid systems. It is shown that the objective function (cumulative oil produced or net present value of the project) is always increased relative to its value in the first generation of the optimization, in some cases by 30% or more. The optimal well type is found to vary depending on the reservoir model and objective function. The effects of reservoir uncertainty are also included in some of the optimizations. It is shown that the optimal type of well can differ, depending on whether single or multiple realizations of the reservoir geology are considered.

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.

A simplified representation of reservoir heterogeneity is applied to the finite difference modeling of horizontal wells operating under two phase flow conditions. The heterogeneity model, referred to as the s-k* approach, was developed and applied in earlier studies for single phase flow scenarios. The model represents heterogeneous permeability fields in terms of a constant background permeability k* and a variable skin s, which captures the effects of near-well heterogeneity. This simplified permeability representation is applied to the modeling of multiple realizations of heterogeneous permeability fields characterized by three different sets of permeability statistics. Comparisons to simulation results using the detailed permeability field and to results using a homogeneous permeability field are presented. The s-k* permeability representation is shown to provide an accurate means of representing reservoir heterogeneity in models of horizontal well performance. Integrated quantities, such as well productivity index, are captured with better accuracy than are more local quantities such as water cut. The s-k* models are more accurate than the homogeneous models in all cases considered. The use of the new modeling approach for the coarse scale modeling of horizontal well performance is also demonstrated.

An efficient semi-analytical solution method appropriate for the approximate modeling of the performance of non-conventional wells operating under primary production is described. The approach, referred to as the s-k* method, accounts approximately for the effects of reservoir heterogeneity and for pressure losses within the wellbore. The accuracy of the method is demonstrated through comparison to detailed finite difference results. Extensive calculations for the performance of non-conventional wells in multiple realizations of log-normally distributed permeability fields are then presented. The effect of wellbore hydraulics on well performance is assessed for several different well geometries and for varying degrees of reservoir heterogeneity. Comparisons between ensemble-averaged quantities (i.e., results averaged over many realizations) and results for a single "average" (homogeneous) reservoir representation are also made. As expected, these calculations differ in general, though for some quantities (e.g., productivity ratio) the homogeneous results give very accurate estimates. The overall approach provides an efficient means for estimating the combined effects of reservoir heterogeneity and wellbore pressure losses on the productivity of non-conventional wells.

A new mechanistic model for gas-liquid two-phase flow with mass transfer through the pipe wall, which accounts for the effects of wall inflow or outflow on wall friction, acceleration and flow pattern transition, has been developed. Four individual flow patterns, including stratified flow, annular-mist flow, bubble flow and intermittent flow, are considered in the mechanistic model. Models for individual flow patterns have been developed. Transition mechanisms among different flow patterns are discussed and new transition criteria are proposed. Predictions by the new mechanistic model are compared with experimental data obtained on a large industrial facility.

A semi-analytical method is presented for the approximate modeling of the productivity of non-conventional wells in heterogeneous reservoirs. The approach is based on Green's functions and represents an extension of a previous model applicable for homogeneous systems. The new method, referred to as the s-k* approach, models permeability heterogeneity in terms of an effective skin that varies along the well trajectory and a constant background permeability k*. The skin is computed through local, weighted integrations of the permeability in the near-well region and is then incorporated into the semi-analytical solution method. The overall method, which can also model effects due to wellbore hydraulics, is quite efficient in comparison to detailed finite difference calculations. Results for the performance of non-conventional wells in three dimensional heterogeneous reservoirs are computed using the s-k* approach and compared to finite difference calculations resolved on the geostatistical fine grid. The new method is shown to provide an accurate estimate of wellbore pressure and production rate, as a function of position along the wellbore, for various well configurations and heterogeneous permeability fields. The possible use of the overall approach in a simulation while drilling (SWD) tool, in which the well path and trajectory are "optimized" using real-time data, is also discussed.

A well model is required to relate the well rate to the well pressure and the block pressure containing the well for the modeling wells in 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 non-conventional wells, when the default procedure available in commercial simulators is not adequate.In this paper we describe a new 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 three-dimension solution for hydraulically fractured wells is obtained by applying Fourier analysis to a two-dimensional solution, therefore this three-dimension solution is easy to obtain when the two-dimensional 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.

To determine the range of validity of pseudo functions for upscaling of multiphase flow with gravity and capillary pressure, the performance of different kinds of pseudo functions is evaluated under different gravity numbers, capillary numbers and upscaling levels. The performance of pseudo functions can be measured as the percentage error between the results of fine grid simulation and the results of coarse grid simulation using pseudo functions. Furthermore, the performance can be shown as the bar graphs of percentage error for varying gravity number, capillary number and upscaling levels. A two-dimensional homogeneous water flooding case and a three-dimensional heterogeneous case have been studied using pseudo functions (Jacks et al.2, Kyte and Berry3, Pore Volume Weighted7, Flux Weighted Potential10, Modified Stone*4, Streamline5), and pseudo function performance has been evaluated under different conditions. Some guidelines for the use of pseudo functions are provided. The results show that (1) all of the pseudo functions have similar performance for the homogeneous case considered; (2) only the Jacks et al. pseudo function can be used in all cases; (3) using pseudo functions (except Jacks et al.) is unreliable, especially for the highly heterogeneous cases.

Fully implicit coupling of surface facilities with a reservoir model can lead to a significant increase in computing time. The focus of this paper is on methods designed to reduce the computational effort. An adaptive implicit method is coupled with a novel preconditioning method to produce an optimal strategy for fully coupled surface facility reservoir simulation.

Reservoir grid generation can be simplified by decomposing the domain into a set of subdomains such that modular grid generation can be performed within each subdomain. The question of how to decompose the domain is addressed. Streamline-potential coordinates form an integral part of our domain decomposition strategy.

Abstract A wellbore/reservoir flow coupling model is developed. It can be applied to a parallelepiped drainage domain consisting of any number of wells with any number of laterals of arbitrary configurations. Reservoir inflow for entire time period (from early transient to late pseudo-steady state) is handled in the coupling model by integrating the instantaneous point source/sink solution over both time and space. Pressure drop along any lateral in the well system includes those caused by wall friction, gravity, and acceleration due to wall inflow or outflow. Effects of wall mass transfer on wall friction and kinetic energy change are taken into account. A computer package has been developed and used to perform a series of sensitivity studies which lead to several important and interesting observations. The coupling model may be applied to determine well productivity, well index, wellbore pressure profile, wellbore inflow or outflow distribution for any time in the well's production/injection life. P. 79

Abstract The fully implicit coupling of surface facilities with a reservoir model can lead to a significant increase in computing time. The focus of this paper is on the development of methods that can minimize the impact of the facility model on the global solution procedure. A local grid refinement preconditioner is developed to accelerate convergence of the full flow field implicit simulation. P. 181

Abstract A brief overview of the 5-year Stanford horizontal wellbore experimental program is presented. Different existing mechanistic flow models, simplified flow models as well as empirical correlations are applied to analyze the 1995, 1996 and 1997 Stanford horizontal wellbore experiments. Comparisons between measurement and prediction are performed, which lead to the conclusion that none of the existing models can predict, with sufficient accuracy, pressure drops for gas-liquid wellbore flow with radial influx.

Abstract Comprehensive, three-dimensional, anisotropic and transient, reservoir/wellbore coupling models are developed for infinite-conductivity and infinite-conductivity wellbores to calculate the productivity of a horizontal well placed in a box-shaped drainage volume. The finite-conductivity well model considers frictional and accelerational pressure drops in the wellbore. It also takes into account the effects of fluid inflow into the well. These models are semianalytical in nature and can be used under transient or pseudosteady state flow conditions. The infinite-conductivity well model is compared with the uniform flux model to show the differences in well productivity calculations. The effects of wellbore pressure drop on well productivity are shown using the finite-conductivity well model. These models can provide reference solutions to calculate well indices or to check the accuracy of a numerical simulator. Because of their analytical nature they are easier and faster to use than numerical simulators. Using superposition in space and time, these models can be run under a variety of constraints. We believe that our finite-conductivity well model is more general than any published previously. P. 191

Abstract Water and gas cresting in horizontal wells are important phenomena in reservoirs that have an aquifer and/or a gas-cap. In practical situations, many reservoirs are produced under supercritical rates and a breakthrough of the displacing phase becomes inevitable. At the beginning of a reservoir simulation study, it is desirable to make an estimate of the breakthrough time and the post-breakthrough behavior. Grid sensitivity runs are also required to obtain the appropriate grid block sizes. An accurate representation of cresting behavior requires a very fine grid, which is not always practical. In this work a procedure was developed to obtain accurate breakthrough times using just coarse grid simulations. The flow equations were written in dimensionless form and important parameters affecting multiphase flow were identified. Simple correlations for a quick estimate of breakthrough time, maximum oil rate and post-breakthrough behavior were derived based on an appropriate set of dimensionless variables and an extensive number of simulation runs. The effects of gridblock size and grid pattern were investigated in detail. Effects of rate, mobility ratio, well drainage area, well height, and end-points and shapes of relative permeability curves were also included. A procedure to derive pseudofunctions either using numerical correlations or coarse grid simulations is also presented. These pseudofunctions can be used to improve the performance of coarse grid simulations. An optimum grid pattern to start a reservoir simulation study is proposed. Application to a real field example shows that the correlations and procedures derived are reliable and accurate, and can be used for quick estimates before starting a reservoir simulation study. Introduction Petroleum reservoirs often have a gas cap and/or an aquifer. In these situations they are subjected to rapid gas and/or water movement towards the well, created by a sharp pressure gradient in the well direction. As the production begins, the interface between the fluids, that is, gas-oil contact or water-oil contact, deforms from its initial plane shape to a cone or a crest. When a field is developed by vertical wells, the shape of the deformation is called a cone, but when it is developed by horizontal wells, this deformation is better described as a crest. The literature reports analytical and semi-analytical solutions for vertical and horizontal wells for critical rate under the steady-state condition reached by the cone or crest with a constant potential on the lateral boundary. There are also analytical and semi-analytical solutions available for infinite acting and closed-boundary reservoirs. Some numerical correlations are also available for predicting breakthrough time and post-breakthrough behavior for supercritical rates. Dimensionless Equations Correlations are more general and reliable when based on the fundamental flow equations, so that the important physical mechanisms are automatically incorporated in the correlations. Here, we consider the case of water cresting. The gas cresting case is treated similarly and is discussed in a later Section. Assuming a homogeneous reservoir with constant viscosity, and neglecting capillary effects, equations for two component, twophase flow, under constant total production rate constraint, can be written as: Oil: (1)

Abstract A novel procedure is proposed for a quick and accurate evaluation of bottomhole wellbore pressure and productivity index. The procedure has led to simple correlations for both horizontal and vertical wells which can be applied to calculate the bottomhole wellbore pressure and productivity index for the entire producing period. The correlations proposed here provide an excellent approximation of the semi-analytical solution at a much reduced computational time. Successful validation of the equations for different well locations and different reservoir geometries has also been achieved. A simple relationship is also presented to estimate the time to reach pseudo-steady state. In addition, change with time of the horizontal well shape factor is shown. The stabilized shape factor for a horizontal well is decomposed into different components which can be easily evaluated by the equations provided. Values for the stabilized well shape factor for various reservoir dimension ratios of practical interest have also been computed and tabulated. Introduction Bottomhole wellbore pressure and productivity index are two important quantities for both horizontal and vertical wells. Approaches reported in the literature to determine these two quantities for horizontal wells can be divided into the following two basic types: Integral transformation (IT) (such as Laplace Fourier transformation) and separation of variables (SOV) methods. Typical examples for this type of approach include Goode & Thambynayagam, Kuchuk et al., Carvalho & Rosa, Aguilera & Ng, Liu & Yang, Ouyang.–Instantaneous point source/sink superposition and the Newman's product method. This kind of approach has been used by a number of researchers, including Gringarten & Ramey, Clonts & Ramey, Daviau et al., Babu & Odeh, Rosa & Carvalho, Ozkan et al., Ouyang, and Economides et al. Although both approaches have been successfully applied to solve a lot of practical reservoir engineering problems, they take nontrivial (and usually long) computing time to get results. As a consequence they are impractical for routine reservoir engineering applications, such as well performance modeling. A novel procedure is proposed here to develop simple correlations for bottomhole wellbore pressure for both horizontal and vertical wells. Problem Formulation Consider a horizontal or a vertical well in a homogeneous but anisotropic parallelepiped reservoir. The reservoir has dimensions of xe, ye, ze in x, y, and z directions, respectively. Its outer boundaries can be either constant pressure or impermeable (Fig. 1). Assume that the formation properties are independent of pressure and reservoir fluid is slightly compressible and of constant compressibility. With these assumptions the partial different equation describing isothermal fluid flow in the reservoir can be expressed as: (1) Coordinate transformation, such as Besson, can be easily applied to transform an anisotropic reservoir to an isotropic one. Hence, only an isotropic reservoir will be considered in the present paper unless otherwise specified. The results presented in the paper can also be used for anisotropic reservoirs by means of the afore-mentioned coordinate transformation. For an isotropic porous media, Eq. 1 becomes (2) P. 535^

The commonly used Cartesian and hybrid local grid refinements have the disadvantage that the base grid is always Cartesian, which makes it difficult to align along reservoir heterogeneities. Also, the locally orthogonal Voronoi and A;-orthogonal perpendicular bisector (PEBI) grids are only useful for isotropic or limited anisotropic permeability distributions. The finite-difference approach proposed in this paper overcomes limitations of existing flexible gridding schemes in modeling full, anisotropic and asymmetric permeability tensors and permeability heterogeneity. The new scheme assumes uniform properties inside control-volumes and it can be used for control volumes formed around vertices of triangles in two-dimensions and tetrahedrain three-dimensions. It can thus be used with Voronoi grids in two-dimensions, and median (i.e. CVFE) and boundary adapting grids (BAG) in two- and three-dimensions. Applications of Voronoi grids in three-dimensions are limited due to geometrical considerations. The proposed method can also be used with grids constructed to align along major reservoir heterogeneities, wells and streamlines. Several applications of this method for two- and three-dimensional problems are presented.

Abstract A novel procedure is proposed for a quick and accurate evaluation of bottomhole wellbore pressure. The procedure has led to simple correlations for both horizontal and vertical wells which can be applied to calculate the bottomhole wellbore pressure and productivity index for the entire producing period. The correlations proposed here provide an excellent approximation of the semi-analytical solution at a much reduced computational time. Successful validation of the equations for different well locations and different reservoir geometries has also been achieved. Approach Bottomhole wellbore pressure and productivity index are two important quantities for both horizontal and vertical wells. For a horizontal or vertical well in a homogeneous but anisotropic reservoir, analytical equations can be obtained for bottomhole wellbore pressure drawdown once pseudo-steady state is reached. Practically, simple equations to evaluate the bottomhole wellbore pressure drawdown in early and transient time are as important as those for late time when the pseudo-steady state flow occurs. Although either the integral transformations (IT) (such as Laplace and Fourier) and separation of variables (SOV) method or the instantaneous point source/sink superposition and the Newman's product method may be used to determine the wellbore pressure drawdown, they take nontrivial (and usually long) computational time to get results, as a consequence they are impractical for routine reservoir engineering applications, such as well performance modeling. Therefore, a novel procedure is proposed here to develop simple yet accurate correlations to evaluate bottomhole wellbore pressure and well productivity index. The idea is very simple and the procedure consists of the following three steps: Step 1. Use the instantaneous point source/sink solution (IPSS) approach to compute wellbore pressure drawdown at different times; Step 2. Design appropriate functional forms which reflect characteristics of the wellbore pressure variation with time; Step 3. Apply a nonlinear regression method to determine the relevant correlation coefficients. Results Correlations for Horizontal Wells. Two functional forms of the wellbore pressure drawdown of horizontal wells, which are valid for the entire time period (from early time to transient and then to late-time long after pseudo-steady state is reached) have been tried for correlating the wellbore pressure drawdown data for different horizontal wells: Form I: (1) Form II: (2) where (3) The values of correlation coefficients, a, b, A and B, for different horizontal wellbore lengths, well locations and reservoir dimension ratios are given in Ref. 2. It was found that both functional forms can provide accurate predictions of the wellbore pressure for the entire time range. P. 287^

Abstract This paper presents a detailed investigation on the reliability of the dynamic pseudo functions used to up-scale flow properties in reservoir simulation. A theoretical study and a real field application are used to evaluate Kyte and Berry (1975). Stone (1991), and a new flux weighted potential (FWP) method. A derivation of Stone's method is presented and shown to use an inconsistent set of equations. Stone's analytical example was used to illustrate how pseudo relative permeabilities that exhibit non-physical behavior may still give acceptable results, but this success can disappear with changes in boundary conditions. The pseudo functions from the field application were not able to match the 2D simulations from which they were calculated, even when a different pseudo function was used for each coarse grid block. Improvements were obtained when directional pseudo functions were used, but still the results were not satisfactory. Similar results were found when comparing fine and coarse grid 3D simulations for a quarter of a five-spot pattern in this field. The results presented in this paper suggest that dynamic pseudo functions, as applied here and as commonly used in industry, may not be an adequate approach to up-scaling. The possibility of large errors and the difficulty in predicting when they may occur make the use of pseudo functions unreliable. Introduction Oil and gas reservoirs are very complex systems in which rock and flow properties vary at all scales (pore to reservoir scale). Rock properties (e.g. porosity and absolute permeability) and flow properties (e.g. relative permeability and capillary pressure) show variations that can be significant to oil recovery at scales below the size of common simulation grid blocks. One of the most important problems in reservoir simulation is that of accurately accounting for such small scale variations. In addition to the low resolution, coarse grid solutions can be strongly affected by numerical diffusion. Several pseudo function techniques have been proposed to reproduce fine grid results (including detailed descriptions of heterogeneity and with minimum numerical diffusion) using the typical coarse grids of field simulations. The dynamic pseudo functions of Jacks et al. (1973) and Kyte and Berry (1975) were designed with the purpose of reducing the dimensionality of a problem (e.g. including 3D effects within a 2D model), and in general, reducing numerical diffusion in coarse gridded models. Simple rules to average absolute permeability were also included in these methods. More recently, up-scaling due to heterogeneities has become the main issue. It was realized that heterogeneities not included in coarse grid simulation could have a major impact on oil recovery predictions. New methods were presented by Hewett and Berhens (1991), Stone (1991), Beier (1992) for dynamic pseudo functions, and by King (1989) and King et al. (1993) for single and two-phase renormalization. The main problem with renormalization is the systematic error that is introduced due to the artificial boundary conditions (e.g. no flow or constant pressure) imposed on each block at each renormalization step. Hewett and Berhens (1993) and Yamada and Hewett (1995) have explained the origin of this error as the forcing of incorrect flow paths or streamlines through the rescaled grid blocks. They have shown that the magnitude of the error can be large and yield incorrect up-scaling (also see Malick and Hewett (1995) for a quantification of the error). The errors in multiphase renormalization are larger than in single-phase since, in addition to incorrect boundary conditions in pressure, multiphase flow renormalization imposes incorrect boundary conditions for fractional flow. The method of Durlofsky et al. (1994) tries to remedy the problems with the renormalization technique. In this method, a non-uniform coarse grid is defined using the flow velocity from a single-phase flow solution (i.e., trying to honor the streamlines), and up-scaling using periodic instead of no flow boundary conditions. Significant errors are still introduced by the artificial boundary conditions (Yamada 1995), and only single-phase flow was up-scaled properly. P. 9

Abstract Fractured reservoirs contain a large fraction of the world supply of oil. For viscous crudes, steam is the most successful technique and field tests indicate that steam has the best potential to recover significant amounts of oil from fractured reservoirs. Unfortunately, there has been little laboratory work done on steam injection in such systems. The experimental system discussed here was designed to understand the mechanisms involved in the transfer of fluids between the matrix rock and the fracture as a result of steam injection. Both continuous and cyclic steam injection experiments were performed on a fractured laboratory system. Saturations were measured in-situ both in the fracture and the consolidated matrix by a CT scanner. The results indicated that there was no steam saturation in the matrix, and that conduction was the dominant heat transfer mechanism. Numerical simulations were used to model both continuous and cyclic steam injection experiments. To model heat losses, heat loss models in the simulator had to be adjusted based on the analysis of the heat losses from the laboratory system with analytical models. After this adjustment, the results from the simulations agreed well with the experiments. When pressure cycling was applied in the simulations with no external heat losses, a considerable amount of steam saturation was observed in the matrix. While the experiments were done with water and steam, simulation runs were also performed for the laboratory system with oil present. Again, steam only flowed in the fracture. Oil recovery was found to be mainly caused by water imbibition into the matrix and heat conduction. Results of this work should be useful in modeling matrix/fracture transfer in dual porosity thermal models. Introduction Fractured reservoirs are estimated to contain 25-30 % of the world supply of oil. Steam injection is required for most of the reservoirs containing heavy oils and tars. There have been quite a few field studies on steam injection for fractured systems (Sahuquet and Ferrier (1982), Britton et al. (1982), Stang and Soni (1987), Closmann and Smith (1983), Duerksen et al. (1984), Couderc et al. (1990) and Hartemink et al. (1995)). Most of these applications were successful or promising. Thomas (1964), Lesser et al. (1966), Abdus Satter (1967) and Satman (1988) presented theoretical models for conduction heating of formations by injecting a fluid through a high permeability streak or fracture. Van Wunnik and Wit (1992) developed a detailed analytical model to study the improvement of gravity drainage by steam injection in a fractured reservoir containing heavy oil. Pooladi-Darvish et al. (1994) developed analytical solutions for heat flow and nonisothermal gravity drainage from a block surrounded by fractures filled with steam. In addition to these analytical models, there have been some numerical simulation studies for steam injection in fractured systems (Pruess and Narasimhan (1985), Lee and Tan (1987), Chen et al. (1987), Pruess and Wu (1989), Nolan et al. (1980), Abad and Hensley (1984), Lin (1988), Briggs (1989), Jensen et al. (1992) and Oballa et al. (1993)). P. 591

SPE Members Abstract The numerical study of water coning requires an accurate spatial discretization around the wells which may lead to a large number of grid blocks. Local refinement is therefore desired since it allows for a coarse underlying grid without loosing numerical resolution around the wells. This paper investigates the applicability of flexible gridding to study coning in both vertical and horizontal wells. These techniques can also be used for multiwell, full field coning studies. A simple sector geometry including a producer and an injector is used. A two step approach is then followed: in the first case, it is assumed that the effects of heterogeneities are negligible and average petrophysical properties are used. In the second part, the heterogeneity of the system is also considered. For both cases, a fine grid simulation on the entire domain is used as the reference response. Various areal hybrid geometries around the producer are then investigated and their response is compared with the reference cases. For the vertical producer, the hybrid geometries are based on a coarse Cartesian grid with either a cylindrical or a Cartesian refinement around the well. For the horizontal producer only an areal Cartesian refinement is considered. In this case though, the hybrid geometry becomes particularly interesting for wells that are not aligned with the overall Cartesian grid. For these cases, the hybrid gridding techniques allow a local Cartesian refinement that is aligned with the well. For all cases investigated, the fine grid results can be re-produced by the hybrid geometries provided that an adequate vertical discretization is used. For both vertical and horizontal producers, the various responses suggest that grid refinement is more important than its geometry. Introduction In standard finite difference simulation the choice of vertical and areal discretization of the domain can have a substantial effect on the numerical results. This is particularly true when steep saturation gradients are present and the underlying grid is too coarse to capture the details of the displacement. A classical example of this problem is the coning phenomenon in wells in which the preferential flow of a phase (usually water or gas) can reduce or even completely inhibit the production of oil through the formation of a 'cone' around the well. In order to include coning in large grid simulations, the usual approach has been to develop either coning correlations or pseudo-functions. In the later case, single well fine grid simulations are performed to back out the pseudo relative permeability and capillary pressure functions so that the coarse grid simulation will match the observed well behavior. For the special case of a homogeneous reservoir some analytical work and practical formulas are also available. These techniques lack generality and are prone to error when applied to multiwell configurations and heterogeneous reservoirs. In recent years, hybrid and refined gridding techniques have been developed into powerful tools for accurate representation of local behavior of single wells in full field models without having to modify the underlying coarse grid. This idea is pursued in this study by comparing the numerical performance of several hybrid grids and locally refined Cartesian geometries on the demanding coning problem. Problem Formulation and Methodology The following situation is considered: a simple rectangular reservoir with a flat top and a bottom aquifer is produced using two wells, an injector and a producer. The vertical injector is completed in the aquifer, and remains unchanged throughout the study. The producer, on the other hand, is allowed to be either vertical or horizontal. P. 397^

Abstract This study addresses the problem of transferring uncertainty in the geological model through the flow simulation model (the comprehensive simulator - CS) up to the reservoir production forecasts. The uncertainty in the geological model is characterized by the differences between many equally probable reservoir descriptions that honor available data and that are generated by geostatistical stochastic simulation. The uncertainty in the production forecasts can be characterized by the variability or spread in the responses derived from the flow simulation, or more precisely by the probability distribution of response variables, such as cumulative oil production, oil recovery, breakthrough time, cumulative water-oil ratio, etc. A new methodology to perform this transfer of uncertainty is summarized. Each reservoir description is first ranked using a fast simulator (FS) rather than the comprehensive flow simulator (CS). A few selected descriptions are then processed through the CS to generate an approximate probability distribution of the reservoir production response. This new methodology is applied with a coarse grid as the FS model to a two-phase, three-dimensional field scale problem with five producing wells, 8×10 6 STm3 original oil in place, and a very large and active aquifer. The impact of uncertainty in absolute permeability and porosity on the field cumulative oil production is considered over 7.8 years of production. Good results are achieved with a 85% to 93% reduction in the computer time over the use of the CS alone. The favorable results achieved show that extension of the proposed approach to other types of problems is worthwhile. Introduction A reservoir flow simulation study combines data from many different sources and processes the data through a complex nonlinear system of equations to generate the reservoir production forecast required for economic analysis. These data include geological data, seismic data, petrophysical data, well data, production data, etc. A problem with these simulations is that the data are not exact and only a small sample is usually available, thus giving the results inherent uncertainty. Uncertainty is the "lack of assurance about the truth of a statement or about the exact magnitude of an unknown measurement or number". Uncertainty involves the risk of making an incorrect decision because the estimates do not agree with reality. Thus, uncertainty is associated with the economic risk analysis of the reservoir production forecast and is a central concern in the decision making process. Because uncertainity in the production forecast results from the interaction of the uncertainties in all sources of data, as well as from the assumptions in the numerical flow model, a procedure to transfer uncertainty through the numerical flow model up to the production forecasts is highly desired. A general procedure is available, and some approximation have been proposed. However, the available general procedure requires such a large amount of computational work that it is not a viable option for most flow simulation studies. Moreover, the approximations proposed by previous investigators cannot handle all the non-linearities of a flow model or are too arbitrary to be of great value. A new systematic methodology for the transfer of uncertainty, requiring reasonable computational work, is applied on a three-dimensional field scale problem. All computer times reported are for a DEC-5810 with 32MB of memory and 1 processor. P. 47^

Abstract The flow rate at each interface between gridblocks is a function of an "effective" permeability of the connection. The upscaling procedure consists of evaluating average connection permeabilities based on a smaller scale description of the reservoir. While several procedures have been reported in the literature to scale up permeability values for Cartesian grids, these procedures have not been extended to flexible grids. The Voronoi grid allows the specification of each gridpoint independently; therefore, it provides flexibility to represent wells and reservoir heterogeneities. This paper describes an upscaling procedure of permeability values at the interface of Voronoi gridblocks. The description of the heterogeneity, e.g., a stochastic image, is independent of the grid geometry. The power law average of the permeability values at each connection is computed automatically by the computer code. The optimum value of power law coefficient is determined by comparing fine and coarse grid results for selected images and well configurations. Several well configurations were investigated with different grid geometries (Cartesian, hexagonal locally refined hexagonal) rising the same description of the reservoir. The optimum value of the power law coefficient was almost the same for all combinations of grid geometries and well configurations (between -0.5 and 0). Also several stochastically generated images of the reservoir, all honoring the same data, were analyzed. Again the optimum value of the power law coefficient was almost the same for all images. The hexagonal grid yields more robust results because each gridblock has more neighbors than the Cartesian grid. The use of local grid refinement around wells reduces the overall error. Introduction The flow rates at each interface (connections) between gridblocks are functions of an "effective" permeability of the connection. This is a practical engineering approach to solve a problem that is very complicated at the pore scale. Therefore, all procedure" to compute average values for permeability aim at obtaining "reasonable" solutions instead of "exact" solutions. For convenience, permeability values are conventionally specified at each gridblock; however, front a reservoir simulation point of view, they are in fact a property of each connection between gridblocks. In some cases, different connections of the same block have substantially different values of permeability, and any attempt to provide a single value for the entire block will miss this characteristic. Furthermore, the average permeability of a connection is a function of the flow pattern and the physical processes. Hopefully, if the grid is fine enough and the averaging procedure adequate, a single value can be computed for each connection to simulate the flow configurations of interest within acceptable errors. An upscaling procedure consists of evaluating average connection permeabilities based on smaller scale description of the reservoir, which may be, for example, a stochastic image (Journel and Alabert). The upscaling procedures that have been proposed in the petroleum literature can be divided into two types. The procedures of the first type attempt to develop a correlation between the small scale data and the average permeability of the gridblock or connection (Warren and Price, Journal et al. Duetsch, Gomez-Hernandez and Journel, King et al. Once the correlation is calibrated, it can be used to compute block permeabilities for different descriptions of reservoir. P. 291^

SPE Members Abstract Numerical reservoir simulation commonly divides naturally fractured reservoirs into matrix and fracture systems. The high permeability fractures are usually entirely responsible for flow between blocks and flow to the wells. The flow in these fractures is modeled using Darcy's law and its extension to multiphase flow by means of relative permeabilities. The influence and measurement of fracture relative permeability for two-phase flow in fractured porous media have not been studied extensively. Transfer of fluids between the matrix and fractures is known to be one of the most important mechanisms of fluid movement in these reservoirs. Therefore, matrix/fracture fluid transfer, matrix and fracture two-phase flow, and interactions among them must be better understood for accurate simulation of oil recovery from naturally fractured reservoirs. Experimental and numerical work on two-phase flow in fractured porous media has been initiated. The process considered is oil displacement by water in fractured porous media. Fine grid simulations of two-phase flow through a fractured core were performed in order to design the experiments and to study the effects of several key variables. These variables included fracture relative permeability, matrix/fracture capillary pressure, and matrix wettability. This paper presents the results obtained from the sensitivity analysis using fine grid simulations. The numerical computations show that fracture relative permeabilities and wettability become important only at low capillary numbers. A dimensionless capillary number at which the fracture relative permeabilities become an important factor was estimated to be 20. The concept of a limiting capillary number may be extended to determine when fracture relative permeabilities are of influence in field scale simulations of naturally fractured reservoirs. Fracture capillary pressure was observed to have a negative effect on water imbibition, reducing the efficiency of the oil recovery. At high capillary numbers, i.e. when capillary forces are more important, water moves faster inside the matrix than in the fracture which results in high oil recovery. At lower capillary numbers, i.e. when viscous forces are more important, water channels rapidly through the fracture with a low recovery efficiency. Comparison of simulations of similar fractured and unfractured systems shows that at high capillary numbers the recovery is higher for the fractured core. At low capillary numbers the unfractured core exhibits higher recovery. Introduction Dual-porosity, dual-permeability formulations are commonly used to model multiphase flow in naturally fractured reservoirs (NFR). The discretized flow equations for this type of model are as follows: (1) where the subscript denotes a given phase (oil, gas or water), m denotes matrix, and f stands for fracture. P. 605^

Abstract In order to get desired accuracy in reservoir simulation, one must refine the grid in certain regions, for example, the regions around wells. When this is done in conventional simulators, unnecessary refinement occurs in some regions. This problem is avoided by the use of "Local Grid Refinement" (LGR) techniques, but the resulting matrix is large and of irregular sparsity. In this paper we address several practical aspects of some methods of LGR, including practical aspects of some methods of LGR, including the use of hybrid grid, and then present three highly efficient methods of solving resulting problems with some Domain Decomposition (DD) techniques. In the first method, a relaxation in the unknowns is performed. In the second one, we use overlapping performed. In the second one, we use overlapping boundaries between the subdomains. In the last one, a coarse grid solution is obtained first and it provides the boundary conditions for the fine grid solution. The LGR and DD methods developed in this work were implemented and tested in a three-dimensional, three-phase, black-oil model. The results are presented in this paper for three examples. It is also presented in this paper for three examples. It is also shown that the DD technique with overlapping boundaries is extremely efficient for solving problems containing one or more LGR regions. The method problems containing one or more LGR regions. The method proposed in this paper can easily take advantage of proposed in this paper can easily take advantage of parallel processors. Also, it is easy to implement the parallel processors. Also, it is easy to implement the proposed techniques in existing simulators. proposed techniques in existing simulators. The LGR and DD techniques discussed in this paper improve the accuracy of reservoir simulation paper improve the accuracy of reservoir simulation around wells and in other regions of high activity. Using these techniques, one not only obtains the correct well pressure, but also the correct saturation near the well. A consequence of this is that WOR and GOR are predicted accurately without well pseudofunctions. pseudofunctions Introduction Reservoir performance forecasting is accomplished by reservoir simulators that solve coupled non-linear partial differential equations describing multicomponent, multiphase flow in the reservoir. In most cases the flow equations are reduced to a system of non-linear algebraic equations by finite-difference techniques and then linearized by the Newton-Raphson method. The resulting linear system has a sparse matrix that can be solved by direct methods for small problems and by iterative methods for larger problems. In the practical application of this tool the number of equations can be more than 100,000 and can even approach 1,000,000. Special gridding and solution techniques are required to take advantage of the physical characteristics of the problem. Recently, there has been a lot of interest in the application of LGR2-10 in reservoir simulation. This technique gives high resolution in regions where the flow rates are high, such as near the wells, without introducing small blocks in regions of low activity. The structure of the Jacobian resulting from LGR can be exploited by the use of DD techniques. P. 245

SPE Members Abstract In this work, a reservoir simulator is used to study the phenomena of gas percolation and gravity segregation, and their effects on reservoir performances. The sensitivity study of reservoir performance to block and time-step sizes shows that for typical solution-gas-drive reservoirs, large en-or's in average reservoir pressure may result from improper control of time-step in the simulator. The comparison of simulation results with the Tarner's method shows that the latter predicts faster reservoir pressure and oil saturation decline, and thus, underpredicts the reservoir producing life and recovery. We propose a new material balance method for predicting the performance of thick, homogeneous, depletion-drive reservoirs. This method accounts for the vertical pressure and saturation gradients, and the secondary gas cap. The thickness of the secondary gas cap can be estimated with good accuracy using an idealized saturation profile. An iterative procedure relates average reservoir pre well pressure. The pressure and saturation at the well are then use to calculate the producing gas-oil ratio. The idealized saturation profile is also used to develop pseudo-functions to simplify the simu-lation of solution-gas-drive reservoir with gravity segregation. Introduction While it is usually possible to do a detailed (and expensive!) three-dimensional, multiblock reservoir simulation study to make predictions for solution-gas-drive reservoirs, it is often instructive to first do some simple material balance calculations Lo make "ball-park" projections. This is especially important because simulation of solution-gas-drive reservoirs involving gas percolation and gravity segregation is a numerically difficult problem that can gobble up large amounts of computer time. The evolution of solution gas and its rapid movement to the top of the reservoir is both the source of numerical difficulties and the justification for assuming that gas movement is essentially instantaneous in performance prediction cal-culations using material balance. The performance prediction of a hydrocarbon reservoir under different drive mechanisms is of interest to any practicing reservoir engineer. A number of studies of solution-gas- drive reservoirs under different conditions have been published. Different conditions that arise during the exploitation of such reservoirs are:Internal gas drive mechanism. Volumetric, undersaturated reservoirs produce by reservoir fluid expansion. The production is attributed to liquid expansion and to rock compressibility, as the reservoir pressure drops down to the bubble point pressure. As the reservoir pressure declines further, oil phase contracts because of the release of solution gas, and the production is due to gas expansion. As gas saturation reaches the critical value, free gas begins to flow, resulting in high gas-oil ratios and low oil recoveries.External gas drive mechanism. In many instances, reservoir pressure in solution-gas-drive reservoirs is maintained by gas injection, and oil is displaced by injected gas. This is referred to as external gas drive mechanism.Gravity segregation (or gravity drainage). For high relief reservoirs with good along-dip permeability, favorable conditions exist for gravity segregation of injected gas or gas released from solution. Gravity segregation is an important factor in attaining high oil recovery from solution-gas-drive reservoirs. Oil recoveries of the order of 60 to 80% can be obtained with effective gravity segregation. Reservoirs with significant gravity segregation also show low producing gas-oil ratios in structurally lower wells. Craft and Hawkins describe these drive mechanisms in more detail. Tarner and Muskat proposed methods to predict the performance of depletion (solution-gas)-drive reservoirs under internal gas drive mechanism, using rock and fluid properties. The assumptions of both methods include negligible gravity segregation forces. Thus, these authors considered only thin, horizontal reservoirs. Both methods use the material balance principle (static) and a producing gas-oil ratio equation (dynamic) to predict reservoir performance at pressures, where gas saturation exceeds the critical value. A more detailed description of both methods appears in Craft and Hawkins The points to note about these prediction techniques are:Time is not a factor in these methods because neither water influx nor gravity segregation are considered. Time history must be inferred from the reserves and well production rates. P. 259^

Abstract Multiphase equilibrium calculations for binary, ternary, and quaternary hydrocarbon-water systems at high temperature were performed using the Schmidtz-Wenzel equation of state. The solubility of water in the hydrocarbon-itch liquid phase and vapor phases is modelled using a constant binary interaction parameter between water and hydrocarbon, while the solubility of hydrocarbons in the aqueous phase is calculated using a temperature dependent binary interaction parameter. Two and three phase equilibrium calculations were performed using the method of successive substitution. A stability analysis using the tangent plane criterion was used to determine the correct number of phases present. At high temperature the solubility of water in hydrocarbon liquids can be quite large. Using the procedure outlined above the solubility of water in the hydrocarbon-itch liquid phase was calculated accurately both for mixtures of water and pure hydrocarbon components, and for mixtures of water and petroleum fractions. The binary interaction parameters used in the hydrocarbon-rich liquid phase and the vapor phase were found to be dependent on the type of hydrocarbon. The interaction parameters for components in the same homologous group such as alkanes, aromatics, and alkenes were found to be almost the same. The solubility of hydrocarbons in the aqueous phase was calculated reasonably accurately using temperature dependent interaction coefficients in the aqueous phase. Above 200 the binary interaction parameters in the aqueous phase were linear functions of temperature. Binary interaction parameters from two phase binary data were found to be quite adequate to calculate three phase multicomponent equilibrium. Introduction Three phase equilibrium calculations with hydrocarbon and water phases are required in the modelling of many reservoir processes. At low temperatures, the solubility of water in the hydrocarbon-rich liquid phase and the solubility of hydrocarbon components in the aqueous phase is small. Thus the hydrocarbon phases and the water phase can be treated as being completely immiscible for the purposes of modelling without introducing significant error. At high temperatures, however, the solubility of water in the hydrocarbon-rich liquid phase can be quite large and the solubility of light hydrocarbons and acid gases in the aqueous phase is appreciable. Hoot and Brady et al. have given experimental data for the solubility of water in pure hydrocarbon components and petroleum fractions at the three phase pressure as a function of temperature. The solubility of water in hydrocarbon liquids increases exponentially with temperature, and at temperatures above 500F, it can be as high as 40 mole %. The solubility of water in hydrocarbons is primarily dependent on the temperature and the paraffin/aromatic characteristics of the hydrocarbons with only a slight effect of carbon number and molecular weights. The solubility of water at a given temperature is the lowest in alkanes and napthenes and the highest in aromatics, particularly benzene (see Figs. 1 and 2). The solubility of water in alkenes lies between its solubility in alkanes and aromatics. The solubility of hydrocarbons in water is considerably less than the solubility of water in hydrocarbons. At high temperature, the solubility of light hydrocarbon components such as methane and acid gas components such as CO2 and H2S can be greater than 10 mole %. The solubility of hydrocarbons in water drops off sharply with increased molecular weight. Given approximate equivalence of molecular weight, alkanes, then alkenes, then napthenes, and finally aromatics are progressively more soluble in water. Several papers in the literature deal with the modelling of hydrocarbon-water systems using equations of state (EOS). Heidemann used the Wilson modification of the Redlich-Kwong EOS to perform three phase flash calculations for hydrocarbon-water systems at temperatures below 300F. The calculation DID procedure was based on the minimization of the Gibbs free energy of the system by integration along the path of steepest descent. Heidemann's scheme used the component material balances and iterated on the component compositions in each phase so that the Gibbs energy of the system was minimized. He tested his procedure by modelling a three component and a six component system. Heidemann used a constant interaction parameter in all phases and was able to match the water content of the hydrocarbon liquid and vapor phases reasonably well. However, the predicted solubility of hydrocarbons in the aqueous phase was orders of magnitude too low. The efficiency of the method was poor, and a large number of iterations were necessary to obtain convergence. P. 733^

Abstract Published imbibition experiments of an advancing fracture water level surrounding a single matrix block are simulated using various double porosity models. These double porosity formulations are inherently unable to represent transverse imbibition with an advancing water level in the fracture as they are based essentially on one dimensional flow. Gross matches of the experimental and double porosity simulations required shape factors that were unrepresentative of the matrix block size. Increasing the fracture water injection rate required increasing the shape factor in order to achieve gross matches, hence a constant shape factor does not give the best matches with experimental data at different rates. A new double porosity transfer function for transverse imbibition is presented based on a nonlinear diffusion equation model for imbibition and the effective exposure time of elements of the matrix block to fracture water. This formulation uses the physical dimensions of the matrix block directly, i.e. no shape factor, and was effective in matching the experimental data without tuning any parameters. A comparison of a cubical single matrix undergoing fracture waterflooding was made between the various double porosity formulations and a fine grid, single porosity simulation. The conventional double porosity simulations did not this fine grid simulation as well as the proposed double porosity formulation. Introduction Naturally fractured reservoirs present many challenges from a numerical modeling point of view. The bulk of the storage of a naturally fractured reservoir resides in the matrix blocks while the fracture network makes up the dominant reservoir flow paths. The overall performance of any naturally fractured reservoir should then be quite sensitive to rate of fluid exchange between the storage body and the reservoir flow paths. This exchange, the matrix/fracture transfer function, is a critical component of any mathematical model used for the simulation of these reservoirs. In this paper, we will present a model for incorporating the effect of transverse imbibition into the matrix/fracture transfer function. We will test the proposed model with two cases, the published experimental work of Kleppe and Morse and a fine grid simulation of a cubical matrix block undergoing fracture waterflooding. Comparisons of the double porosity imbibition modeling methods presented by Litvak and by Thomas et al. will also be made to these two test cases. Warren and Root presented equations for unsteady-state single phase flow in double porosity reservoirs. Their double porosity domain assumes a continuous. uniform fracture network oriented parallel to the principal axes of permeability. The matrix blocks in this system occupy the same physical space as the fracture network and are assumed to be identical rectangular parallelepipeds with no communication between matrix blocks. Matrix blocks are also assumed to be isotropic and homogeneous. Assuming pseudosteady-state flow between the matrix elements and the fracture network, Warren and Root presented analytical solutions of these double porosity equations for pressure build-up tests. Their matrix/fracture transfer rate is controlled by a geometrical shape factor which is a function of the surface-volume ratio of the matrix blocks. Extension and numerical solution of the Warren and Root double porosity model to multiphase flow in three dimensions was done by Kazemi et al. Attempts to match water imbibition into artificially fractured cores by Kazemi and Merrill showed some success if the matrix block was divided into subdomains. As capillary pressure curves were not measured on the cores, capillary pressure was a parameter that was modified to help fit the recovery curves. cater enhancements of the matrix/fracture transfer function by Gilman and Kazemi included a more realistic gravity potential between the matrix and fracture and the ability to account for fluid transfer to the matrix system due to an imposed pressure gradient in the fracture. The improved gravity transfer was accomplished by gridding the matrix block into subdomains. P. 155^

SPE Members Abstract Phase injectivities for multiphase injection processes are studied using an isothermal black-oil numerical model that properly treats wellbore/reservoir interaction. A fully implicit technique for the coupling of wellbore and reservoir flow equations is described. The effect of gravity segregation in the wellbore is taken into account to simulate cocurrent water and gas injection. It is shown that the current techniques used in reservoir simulation for assigning phase injectivities yield inconsistent results when wellbore phase redistribution takes place. The effect of multiphase flow in the well and in the reservoir during well testing is also investigated. The fully coupled wellbore/reservoir model is employed to study the effect of wellbore phase segregation on buildup pressure response. pressure response. Introduction Multiphase injection of fluids in reservoirs may occur in a variety of oil field operations such as steam and cocurrent water and gas injection processes. In enhanced recovery processes, special techniques are usually required to control and modify injection profiles. Laboratory experiments reported by Elson have showed the effect of phase separation in the wellbore. Uneven injection profiles in the perforated intervals have been clearly observed. Field observations have also shown this behavior with the lighter and less viscous phase being preferentially injected into the top of the formation while the heavier phase goes into the lower part. phase goes into the lower part. In reservoir simulation, phase injectivity is commonly handled by means of simplified approaches which usually ignores the gravity segregation in the wellbore. Multiphase flow in the wellbore and in the surrounding formation may also occur during well tests. Anomalous pressure buildup behavior, as a result of wellbore phase pressure buildup behavior, as a result of wellbore phase segregation, has been observed in high GOR wells. A general review of well test analysis with multiphase flow was presented by Raghavan. However, the effect of wellbore phase segregation was not taken into account. An attempt to account for this effect in pressure buildup analysis was made by Fair. In his work, an empirical relationship for pressure change was introduced in the wellbore storage equation in order to describe phase redistribution in the well. Recently, Winterfeld presented a model that treats rigorously multiphase flow in the wellbore during pressure buildup. The model solves simultaneously for transient multiphase flow in the wellbore and in the reservoir. Counter-current two-phase flow in the wellbore during shut-in is allowed by the use of semi-empirical expressions to account for phase to phase viscous forces. To allow for phase segregation in the simulation of multiphase injection processes and multiphase well tests, a wellbore flow model must be coupled with the reservoir model. Several papers dealing with the simulation of transient multiphase flow in pipes have been published in the technical literature. Liles and Reed developed a semi-implicit drift-flux model for simulating unsteady state two-phase flow in pipes. Later, Millers presented a similar model to study wellbore storage effects during geothermal well testing. Sharma et al. simulated transient gas-oil flow in pipelines using a black-oil type model. A thermal pipelines using a black-oil type model. A thermal compositional model for simulating transient gas-liquid flow in natural gas pipelines was reported by Kohda et al.. P. 123

SPE Members Abstract Published imbibition experiments of an advancing fracture water level surrounding a single matrix block are simulated using a fine grid single porosity model and various double porosity models. The fine grid simulations show that a stationary water saturation profile quickly develops and advances in the matrix at the same rate as the fracture water level. A new double porosity transfer function for imbibition dominated matrix/fracture fluid exchange is presented based on stationary profile solutions of the fractional flow equation. This transfer function models the imbibition recovery for these advancing water level experiments better than the conventional double porosity imbibition formulations. It is shown that the stationary profile transfer function is best suited to systems where the time to develop the stationary profile is short relative to the length of the waterflood. The experimental data was also simulated using a diffusion equation with a nonlinear diffusion coefficient in combination with a moving boundary condition as the imbibition model. A constant diffusion coefficient based upon water relative permeability and capillary pressure gradient values at 1 - Sor matched the experimental results as well as the nonlinear diffusion coefficient, but required much less computer time. From analyzing two different diffusion type equations as imbibition models, we show that countercurrent imbibition is not a likely recovery mechanism for this type of advancing water level imbibition experiment. Introduction Naturally fractured reservoirs present one of the most difficult problems in reservoir simulation. The fracture network may be spatially difficult to define while the transfer of fluids from the matrix blocks to the fracture network may be qualitatively described by many different physical mechanisms. This paper will discuss the exchange of fluids via capillary imbibition in a matrix domain of regular, dispersed parallelepipeds within a continuous fracture domain. We will present fine grid simulation of the experimental work of Kleppe and Morsel and present a transfer function that models their single block waterflood better than existing double porosity formulations. The double porosity approach formulated by Barenblatt et al. for single phase flow considered the naturally fractured reservoir to be composed of two superimposed continua, a continuous fracture system and a discontinuous system of matrix blocks. The fracture system has a low storativity and high conductivity while the majority of the oil resides in a matrix block system of low conductivity. The transfer of fluids between these two systems, as described by Barenblatt, assumes pseudosteady-state flow between the matrix blocks and the fracture system. The matrix/fracture transfer function can then be expressed as Darcy's Law with an appropriate geometrical factor that describes the characteristic length and flow area between the matrix blocks and the fracture system. Warren and Root" presented an analytical solution for single-phase, unsteady-state, radial flow in a naturally fractured reservoir and introduced the dual porosity concept to petroleum engineering. They also derived a formula for the shape factor for parallelepiped matrix blocks within an orthonormal fracture system of one, two or three dimensions. The formulation of Warren and Root was extended to a multiphase system by Kazemi et al., who numerically solved for flow in three dimensions. As with the Barenblatt formulation, Kazemi's transfer function assumes pseudosteady-state flow based on potential differences between the matrix node center and the fracture node center within a gridblock. In a later paper Kazemi and Merrill used their simulator to model water imbibition into artificially fractured cores. Most of their simulations involved gridding the matrix block into more than one computational unit, thereby improving the model's performance through enhanced matrix block definition. With a single matrix block per grid block, they presented simulations at two different water injection rates for cores initially saturated with 100 percent diesel oil and for cores with n-decane replacing diesel oil. P. 509^