A history-dependent model for saturation functions, combined with a three-dimensional, three-phase, semi-implicit reservoir simulator, has been developed. In water-coning simulations with variable rates, for waterflooding in the presence of free gas saturations, and for gas-cap shrinkage, use of hysteresis in saturation functions shows results significantly different from those obtained by conventional methods. To some extent, the model is based upon remembering the saturation history of the reservoir. In doing this, smooth transitions of both relative permeabilities and capillary-pressures from permeabilities and capillary-pressures from drainage-to-imbibition or imbibition-to-drainage states are allowed. In addition, the effect of trapped gas or oil saturations on relative permeabilities and capillary pressures is accounted for. Tests of the model indicate that simulation with hysteresis is a stable-procedure requiring little more computation time and storage than normal simulations. In addition, results of these tests agree qualitatively with experimental and field results. Introduction Present-day reservoir simulators have allowed investigation of complex recovery schemes and production schedules. Although simulators can production schedules. Although simulators can handle such problems numerically, most treat saturation functions in a simplified manner. For example, only one set of saturation functions may be used for initialization and/or simulation in a particular part of the reservoir. It is assumed that particular part of the reservoir. It is assumed that saturation changes occur in a given direction - drainage or imbibition-for most of the simulation. Cutler and Rees pointed out that hysteresis in capillary pressures may affect well coning behavior. Other authors have shown that hysteresis in relative permeabilities is important in the correct prediction of reservoir behavior. Unless treated prediction of reservoir behavior. Unless treated more realistically, the history dependence of saturation functions could cause significant errors in reservoir simulation. This paper describes a reservoir simulation technique in which saturation-function hysteresis is accounted for. A model for hysteresis is incorporated, permitting smooth transitions in either direction between drainage and imbibition relative permeability and capillary pressure curves as observed in experimental data. Including this hysteresis model allows the simulator to predict more realistically many reservoir situations. THE HYSTERESIS MODEL The hysteresis model allows both capillary pressures and relative permeabilities to range pressures and relative permeabilities to range between imbibition and drainage curves via intermediate "scanning" curves. Experimental data are required only for the bounding imbibition and drainage functions since the model provides an interpolative scheme for arriving at the intermediate values. However, regression parameters are incorporated allow a closer fit with experimental scanning states, should these data exist. The model also allows the use of analytical curves for the bounding relative-permeability functions, for which data may not exist. The hysteresis model has been designed so that saturation functions derived from the hysteresis algorithm approach physical reality. To this extent, the existing experimental data have been used as the basis for the model. The following sections describe these data and the associated procedures for calculating hysteretic relative permeabilities and capillary pressures. Further details and equations are given in the Appendix. CAPILLARY HYSTERESIS Capillary hysteresis is characterized by bounding imbibition and drainage curves and intermediate scanning curves, as shown in Fig. 1. SPEJ P. 37

SPE Member Abstract The Ninth SPE Comparative Solution Project presented in the following paper provides a reexamination of black-oil simulation based on a model of moderate size (9,000 cells) and with a high degree of heterogeneity provided by a geostatistically-based permeability field. Nine participants provided data for the comparison which is based on a dipping reservoir with twenty-five somewhat randomly placed producers and a single water injector. Results showed that significant agreement could be achieved for this problem on the basis of total production rates, saturations, and reservoir pressures. On the other hand, rates for some individual wells did show variations of as much as 30% due to differing treatments of well flowing bottomhole pressures. All participants were able to simulate the study in fewer than sixty time-steps with an average of 4-5 Newton iterations per step. In addition, the results showed that this moderate-sized problem could be simulated in only a few minutes in a workstation environment for the two plus years of data. Introduction The SPE Comparative Solution Projects have to this date comprised studies over the past fifteen years which involved varying aspects of reservoir simulation. The first two comparative solution projects focused on black-oil simulation for Cartesian and radial grid geometries. More recent studies have been devoted to specialized simulations such as compositional, dual porosity, thermal, and miscible; or they have looked at treatment of horizontal wells and gridding. The purpose of the current project is two-fold: to provide an update to previous black-oil based projects and to investigate the complications brought about by a the high degree of heterogeneity in a geostatistically-based permeability field. Description of the Reservoir Simulators From the first comparative solution project to the present, the capabilities of black-oil reservoir simulators have increased tremendously. Currently, the linear equation solutions are generally performed by a preconditioned conjugate-gradient-like method such as incomplete LU factorization preconditioned ORTHOMIN. These solvers have predominantly replaced techniques such as SIP and LSOR for the more difficult non-symmetric cases resulting from highly-heterogeneous and fully-implicit solutions. Well treatments have become more sophisticated with almost all simulators providing some form of treatment of implicit bottomhole pressures for rate constrained wells. P. 135

The state of the art of modeling fluid flow in shale reservoirs is dominated by dual porosity models which divide the reservoirs into matrix blocks that significantly contribute to fluid storage and fracture networks which principally control flow capacity. However, recent extensive microscopic studies reveal that there exist massive micro-and nano-pore systems in shale matrices. Because of this, the actual flow mechanisms in shale reservoirs are considerably more complex than can be simulated by the conventional dual porosity models and Darcy's Law. Therefore, a model capturing multiple pore scales and flow can provide a better understanding of complex flow mechanisms occurring in these reservoirs.Through the use of a unique simulator this paper presents a micro-scale multiple-porosity model for fluid flow in shale reservoirs by capturing the dynamics occurring in three separate porosity systems: organic matter (mainly kerogen), inorganic matter, natural fractures. Inorganic and organic portions of shale matrix are treated as sub-blocks with different attributes, such as wettability and pore structures. In the organic matter or kerogen, gas desorption and diffusion are the dominant physics. Since the flow regimes are sensitive to pore size, the effects of nanopores and vugs in kerogen are incorporated into the simulator. The separate inorganic sub-blocks mainly contribute to the ability to better model dynamic water behavior. The multiple porosity model is built upon a unique tool for simulating general multiple porosity systems in which several porosity systems may be tied to each other through arbitrary transfer functions and connectivities. This new model allows us to better understand complex flow mechanisms and in turn is extended into the reservoir scale considering hydraulic fractures through upscaling techniques. Sensitivity studies on the contributions of the different flow mechanisms and kerogen properties give some insight as to their importance. Results also include a comparison of the conventional dual porosity treatment and show that significant differences in fluid distributions and dynamics are obtained with the improved multiple porosity simulation. Finally a case for reservoir-scale model covering organic matter, inorganic matter, natural fractures and hydraulic fractures is presented and will allow operators to better predict ultimate recovery from shale reservoirs.

A typical tight oil reservoir such as the Bakken has matrix pore sizes ranging from 10 nm to 50 nm. At such small scales the confined hydrocarbon phase behavior deviates from bulk measurements due to the effect of capillary pressure. In addition, compaction of pore space can bring about order of magnitude changes for tight oil formation properties during pressure depletion further exacerbating these deviations. Without considering these facts a conventional reservoir simulator will likely not be able to explain the inconsistent produced GOR observed in the field compared to simulated results. The effect of these inaccuracies on ultimate recovery estimation can be devastating to the underlying economics. This paper presents a compositional tight oil simulator that rigorously models pressure dependent nanopore-impacted rock and fluid properties, such as suppression of bubble point pressure, decrease of liquid density, and reduction of oil viscosity as well as their interactions with pore space compaction. The cubic Peng-Robinson equation of state is used for phase behavior calculations. Capillary pressure is evaluated by standard Leverett J-function for porous media. Modifications to the stability test and two-phase split flash calculation algorithms are provided to consider the capillarity effect on vapor-liquid equilibrium. The simulator can capture the pressure-dependent impact of the nanopore structure on rock and fluid properties. As a result, the problem of inconsistent GOR is resolved and the history matching process is greatly facilitated. It is shown that inclusion of these enhanced physics in the simulation will lead to significant improvements in field operation decision-making and greatly enhance the reliability of recovery predictions.