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