Two numerical models are introduced for simulating three‐dimensional, two‐phase fluid flow and heat transport in geothermal reservoirs. The first model is based on a three‐dimensional formulation of the governing equations for geothermal reservoirs. Since the resulting two partial differential equations, posed in terms of fluid pressure and enthalpy, are highly nonlinear and inhomogeneous, they require numerical solution. The three‐dimensional numerical model uses finite difference approximations, with fully implicit Newton‐Raphson treatment of nonlinear terms and a block (vertical slice) successive iterative technique for matrix solution. Newton‐Raphson treatment of nonlinear terms permits the use of large time steps, while the robust iterative matrix method reduces computer execution time and storage for large three‐dimensional problems. An alternative model is derived by partial integration (in the vertical dimension) of the three‐dimensional equations. This second model explicitly assumes vertical equilibrium (gravity segregation) between steam and water and can be applied to reservoirs with good vertical communication. The resulting equations are posed in terms of depth‐averaged pressure and enthalpy and are solved by a two‐dimensional finite difference model that uses a stable sequential solution technique, direct matrix methods, and Newton‐Raphson iteration on accumulation and source terms. The quasi‐three‐dimensional areal model should be used whenever possible, because it significantly reduces computer execution time and storage and it requires less data preparation. The areal model includes effects of an inclined, variable‐thickness reservoir and mass and energy leakage to confining beds. The model works best for thin (<500 m) reservoirs with high permeability. It can also be applied to problems with vertical to horizontal anisotropy when permeability is sufficiently high. Comparisons between finite difference and higher‐order finite, element approximations show some advantage in using finite element techniques for single‐phase problems. In general, for nonlinear two‐phase problems the finite element method requires use of upstream weighting and diagonal lumping of accumulation terms. These lead to lower‐order approximations and tend to obviate any advantage