As for enhanced oil recovery (EOR) methods in heavy oil reservoirs where the existing techniques find their limits, the solvent-based techniques have attracted extensive attentions recently. CO2 is found to be a prevailing solvent due to its superiority of diluting heavy oil as well as its potential of mitigating greenhouse gas emissions, though its solubility is limited in heavy oil. The addition of light hydrocarbon solvent(s) (e.g., C3H8 and n-C4H10) in CO2 stream facilitates its dissolution in heavy oil, resulting in an enhanced viscosity reduction and swelling effect. Furthermore, water exists in the reservoir not only as the formation water (e.g., connate water and bottom water), but also as the injected water during EOR processes (e.g., water-alternating-gas method). Therefore, it is of fundamental and practical importance to quantify the phase behaviour of alkane solvent(s)–CO2–water–heavy oil systems at high pressures and elevated temperatures. A new alpha function in the Peng–Robinson equation of state (PR EOS) has been developed for water component to improve its prediction for water vapour pressure in a wide temperature range. Subsequently, a new temperature-dependent binary interaction parameter (BIP) correlation for the CO2–water pair in the aqueous phase is proposed by matching CO2 solubility in water. The compositions of both water-rich phase and CO2- rich phase can be predicted with a good accuracy by using the newly developed mathematical model. CO2–heavy-n-alkanes systems are commonly taken as a reference of CO2–heavy oil systems. As such, a new BIP correlation for CO2 and heavy-n-alkanes from n-decane (n- C10H22) to n-tetratetracontane (n-C44H90) has been developed. A pressure-associated optimization approach is applied to optimize a BIP value between CO2 and a specific nalkane at an isotherm, resulting in an obvious improvement to predict the phase behaviour of CO2–heavy-n-alkanes system at high pressures. Subsequently, multiphase boundaries and swelling factors are determined for binary systems of CO2–heavy oil, C3H8–heavy oil, n-C4H10–heavy oil, ternary systems of C3H8– n-C4H10–heavy oil, C3H8–CO2–heavy oil, n-C4H10–CO2–heavy oil, and quaternary C3H8–n-C4H10–CO2–heavy oil systems at high pressures and elevated temperatures. Experimentally, a pressure–volume–temperature (PVT) setup is used to conduct the phase behaviour tests of various compositions of alkane solvent(s)–CO2–heavy oil systems. Theoretically, the PR EOS coupled with a previously modified alpha function is applied to quantify phase behaviour of the aforementioned systems by characterizing heavy oil as six pseudocomponents. Compared with an internal tuning process embedded in WinProp Module (version 2011, CMG), a tuning process developed in this work is demonstrated to provide a better performance of matching the experimental measurements. In addition, it is also found that a good initial value of BIP matrix input in the WinProp Module can significantly improve its regression results. Finally, water is added into the C3H8–CO2–heavy oil systems to simultaneously form three phases (i.e., water-rich aqueous (A), oil-rich liquid (L), and solvents-rich vapour (V)) at certain conditions. Not only have the boundaries between AL and ALV phases, compositions of L and V phases, and their phase volumes been experimentally measured, but also can be theoretically predicted with a good accuracy by utilizing the newly developed techniques.