The Fully Implicit Method (FIM) is widely used in numerical reservoir simulation due to its unconditional stability, which allows for arbitrarily large time steps. However, FIM is computationally expensive per time step, especially for large numbers of components and highly detailed reservoir models. IMPEST (Implicit Pressure, Explicit Saturations, Temperature and mole fractions), on the other hand, is computationally inexpensive, but only conditionally stable. For large-scale heterogeneous models, the allowable stable time step of IMPEST may be extremely small. In the Adaptive Implicit Method (AIM), only a subset of the primary variables is treated implicitly. AIM offers a balance between FIM and IMPEST by employing implicit treatment only when and where necessary. The challenge with AIM is to: (1) find robust and sharp stability criteria that can be used to choose an optimal time step size, and (2) devise an efficient switching algorithm to dynamically label variables in gridblocks implicit or explicit. The objective of this thesis is to formulate, implement, and validate an efficient Thermal Adaptive Implicit Method (TAIM) that is capable of solving the nonlinear system of algebraic equations efficiently using the minimum number of implicit variables for a given time step. TAIM includes the stability criteria for the selection of the time step size and a switching algorithm to label the variables as either implicit or explicit. The TAIM stability criteria are obtained using the von Neumann approach. The derivation of the criteria is obtained using a comprehensive linear stability analysis that takes into account the complex physics being modeled, including mass and heat convection, thermal conduction, phase change, and gravity. The criteria were implemented and evaluated using Stanford's General Purpose Research Simulator (GPRS). The stability criteria serve as an input to a switching algorithm that calculates the maximum possible time step size for which a stable solution for explicit saturation, temperature, and compositions in a gridblock is guaranteed. We implemented two switching algorithms in GPRS, namely, percentage-based and variable-based. In the percentage-based approach, the user specifies the percentage break-up between various implicit schemes. Since an arbitrary percentage break-up in difficult problems can lead to extremely small time steps, thereby rendering the TAIM method infeasible, a more optimal solution is for the user to specify a desired time step size, and let the TAIM switching algorithm decide on the implicit/explicit labeling of variables in the gridblocks. We formulated this algorithm and refer to this approach as the variable-based TAIM. We demonstrate that a TAIM-based approach is a promising technique for the simulation of thermal-compositional displacement processes of practical interest. TAIM simulations are not only more accurate than their FIM counterparts due to reduced numerical diffusion, but TAIM also offers an efficient numerical simulation method that uses time steps comparable to those used in FIM, while treating a large subset of the variables explicitly. In this thesis, we derive the thermal-compositional stability criteria that provide the CFL (Courant-Friedrichs-Lewy) numbers to obtain the maximum allowable time steps. The criteria are used in conjunction with a variable-based switching algorithm that labels the variables as implicit/explicit, which can vary both in space and time. We found the stability criteria to be sharp, i.e., when convection is the dominant mechanism, mild violations of the stability limits lead to unstable solutions both in the saturation and temperature...