The production of heavy oil is gradually becoming more viable, as both the price of production decreases and the price of oil (over the long term) increase. The reserves of heavy oil are far greater than the existing reserves of conventional oil; however, it is often not economically desirable to extract. Heavy oil is characterised as a very viscous and dense substance, which flows very slowly under gravity at atmospheric temperature. As such, it is often necessary to use steam to reduce the viscosity of the oil to enhance production. Novel production enhancement techniques can utilise a light hydrocarbon solvent to further increase recovery. The simulation of these reservoirs is not trivial, and during production, they undergo large changes in temperature and pressure depending upon the production strategy used. Furthermore, the use of a light hydrocarbon solvent can lead to highly asymmetric mixtures and the appearance of a second oleic phase. However, an accurate simulation is necessary due to the high cost of the production and the requirement for reliable predictions for new production enhancement strategies. In this work, the thermodynamic models used to represent the fluid and the methods of flash calculation suitable for multiphase mixtures are studied. The final aim is to develop a reservoir simulator incorporating some of these elements. To accurately represent the fluid phases in equilibrium in a heavy oil reservoir, it is necessary to use an equation of state (EoS). There are a huge number of possible equations of state (EoS’s) available. In this work, a number of these models are compared for heavy oil-related binary systems to evaluate which are best suited. The accuracy of each model is compared, and it is found that if the interaction between the water and hydrocarbon phases is not considered to be important, then the commonly used cubic EoS’s with the van der Waals mixing rules are suitable (optimally using two binary interaction parameters). If more complex interactions are considered important to the simulation, then the cubic EoS’s can be used with the Huron–Vidal-type mixing rules (using the non-random two-liquid (NRTL) activity coefficient model). The more complex EoS’s such as cubic plus association (CPA) or PC-SAFT can also be used and are generally as accurate as the Huron–Vidal type. Another issue associated with thermal simulation is that it is necessary to add an energy balance to the system of PDEs used. It is common to use the temperature as an additional primary variable. The resulting isothermal flash specification can lead to some problems during transient simulation. The fluid in some of the grid blocks may become narrow boiling. The enthalpy of a narrow boiling mixture changes rapidly due to very small changes in the temperature or pressure. This can lead to issues with convergence and potential oscillations during Newton steps. A related issue is when the solution requires more phases than components to meet the energy balance equation, and this cannot be solved when using isothermal flash. This can often occur close to steam injection wells (where components other than water can be stripped away). An approach using the energy directly in the flash calculation (e.g. isobaric and isenthalpic flash) is explored in this work and an algorithm presented which takes only slightly more computational resources than the conventionally used isothermal flash. This is demonstrated to be robust for a number of mixtures, and an approach tailored for the flash calculation of mixtures containing water is described. The issue of multiple coexisting equilibrium phases is also examined. The conventional isothermal flash framework is complex and cumbersome when dealing with more than two or three phases. A new method (modified RAND) is presented based on the related chemical equilibrium problem. This method is primarily examined for multiphase isothermal flash, with a robust implementation described. An extension to other state function-based flash specifications is developed using the new proposed framework. Furthermore, the conventional method of solving an EoS at a specified temperature and pressure is abandoned and a method which co-solves the EoS with the equilibrium equations is described (vol-RAND). A different approach to solving state function-based flash specifications is also considered, where the EoS is solved at the given state function variables resulting in a minimisation without nonlinear constraints. This proposed approach has not yet been tested. A thermal, EoS-based reservoir simulator is developed from an isothermal simulator. The energy balance partial differential equation is added, and the temperature is used as an additional primary variable. The EoS’s compared for heavy oil-related fluids are implemented in the simulator. A multiphase flash algorithm using modified RAND and stability analysis skipping is added to the simulator. Finally, a test case is used to demonstrate that some of the more complex EoS models can be used and the multiphase flash algorithm is suitable. With the thermal reservoir simulation tool developed, it is possible to carry out further comparisons and add more complexity in future work.