How to treat well boundary conditions is a difficult issue when formulating and coding a multiphase numerical reservoir simulator. The difficulty arises because the partial differential equation governing multiphase subsurface flow is of a mixed parabolic-hyperbolic type. Maximum changes in primary variables and mass/heat fluxes occur at boundaries or well nodes. Consequently, these well nodes with small volume lend to be singular, leading to computational convergence problems. The conventional method of well treatment in geothermal or oil reservoir simulators is to use a sink/source term approach and to distribute flow rates by a potential or mobility allocation scheme for a multilayered well. However, this traditional method cannot handle a backflow problem, which may occur in a multilayered well in heterogeneous formations. This paper presents a virtual node method to handle a well bore either as a single node or several computational nodes screened and connected to many neighboring nodes for a multilayered well. The well bore can be vertical, inclined, or horizontal, and the well borehole node is treated in the same way as any other nonwell node for flow calculations. The solution at the well is then obtained by solving mass balance equations for the well node. It is shown that the new method provides a natural, physically consistent, and numerically efficient approach to handling well flow problems. Implementation of this new method for a three-phase flow reservoir simulator is discussed, and three examples are provided. ¿ 2000 American Geophysical Union