Abstract In this research study, we conducted a coupled thermal-stress-fluid flow numerical model on the UTF Phase A SAGD project in order to investigate the fundamental geomechanical behaviour involved in the SAGD process, and to gain insight into the reservoir response to temperature and pore pressure changes. The numerical simulation is carried out by using a self-developed coupled finite element model which incorporates our proposed strain-induced permeability model. The obtained simulation results were compared with the measured data. Introduction Thermal recovery processes involve coupling between heat transfer, multiphase flow and stress/deformation, which has become an increasingly important subject in the petroleum field(1). Particularly, the coupling is crucial in problems such as borehole stability, hydraulic fracturing and injection/production induced deformation of the ground surface during the thermal recovery process in heavy oil or oil sand reservoirs. Numerical modelling of the coupled processes is historically carried out in the areas of geomechanics modelling and reservoir simulation. The former is to compute the stress-strain behaviour; therefore, the deformation. The latter is to essentially model the multiphase flow and heat transfer in porous media. Each of these disciplines simplifies the part of the problem that is not of primary interest. These approaches are unacceptable in situations where the coupling is strong and the changes of porosity and permeability cannot be accounted for by rock compressibility alone. Gutierrez and Lewis(2) extend Biot's theory to multiphase fluid flow in deformable porous media. Based on their formulation, they conclude that the coupling between the geomechanics and the multiphase flow occurs simultaneously. Thus, fully coupled system equations of deformations, multiphase flow and heat transfer should be solved simultaneously. Development of such kinds of fully coupled geomechanics-multiphase flow-heat transfer simulators needs tremendous effort, since the existing FEM geomechanics codes and FDM reservoir simulators cannot be used. Settari and Mourits(3) present an approach to couple the stress-strain behaviour to multiphase flow-heat transfer using porosity as a coupling parameter. The geomechanics module and the thermal reservoir simulator are used in a staggered manner. Pore pressure and temperature changes are calculated from the thermal reservoir simulator and transferred to the geomechanics module. The stress and the displacement changes are then calculated in the geomechanics simulation. An iterative algorithm is used to ensure that the porosity calculated from the geomechanics module is the same as that from the thermal reservoir simulator. The staggered technique employed to solve the coupled system equations allows for the use of the existing geomechanics codes in conjunction with a standard reservoir simulator. Currently, most of the commercial coupled geomechanics-multiphase flow-heat transfer simulators are developed in this way. The disadvantage of these kinds of coupled simulators is that the thermal reservoir module, usually developed using the finite difference method (FDM), cannot accommodate the full permeability tensor, since they adopt the standard discretization scheme, such as 5-spot for 2D problems and 7-spot for 3D problems. In this paper, the development of a coupled geomechanics-multiphase flow-heat transfer simulator using the finite element method (FEM) is described with the use of Galerkin's least squares (GLS) technique(4) to stabilize the saturation equation.