Tremendous resources of heavy oil are located in western Canada. Among the many heavy oil recovery methods proposed, thermal recovery methods which employing hot fluid injection or in-situ combustion, have been conventionally utilized to enhance heavy oil recovery. Temperature profiles in heavy oil reservoirs are important factors for making operation and production plans in thermal recovery processes. A significant factor in these processes is to study the heat energy transfer in the oil formations. In previous studies, heat transfer by conduction has been commonly considered as the main mechanism of heat transfer in porous media. However, this assumption is not reasonable for thermal recovery processes with fluid flow such as steam flooding, hot water flooding, steam-assisted gravity drainage (SAGD) and cyclic steam stimulation (CSS). Heat transfer and fluid flow occur simultaneously during thermal recovery processes, and heat transfer by convection cannot be neglected. Fluid flow motivates convective heat transfer, and increase the rate of energy transfer significantly. Heat conduction is dominated by temperature gradient, while heat convection is dominated by pressure gradient. By integrating conduction and convection, temperature and pressure domain can be coupled systematically. In this study, novel heat transfer models, integrating both conduction and convection, have been developed to describe the one-dimensional transient heat transfer coupled with fluid flow. In the models, the properties of the reservoir and the fluid are integrated into two important parameters, i.e., the thermal diffusivity of a reservoir/fluid system and the thermal convection velocity of the fluid. To derive the analytical solutions of these mathematical models, dimensionless variables are defined to reduce the models to the dimensionless form. After that, variable transformation and Laplace transformation are performed to derive the analytical solution in Laplace domain. By using the table of Laplace transformations, the solution in Laplace domain can be converted to dimensionless analytical solution in real time domain. Numerical simulations by COMSOL Multiphysics are conducted to validate the analytical solutions. Subsequently, case studies under both steady and unsteady flow conditions have been conducted. Satisfactory agreements of the results are achieved between analytical solutions and numerical simulation results. To prove the mathematical models could have practical application in the oil and gas industry, results comparison between the analytical solution and CMG simulation is conducted. A numerical simulation model for transient heat transfer in heavy oil reservoirs during the SAGD process was used for comparison. It is found that the shapes of temperature distributions and propagations of the analytical solution and the CMG simulation have the similar trends. The studies showed good agreement between the test results and those from the CMG simulation. The newly developed analytical solutions provide theoretical guidance for temperature transient analysis (TTA) and fluid injection strategies. These analytical solutions can be used to predict temperature profiles in heavy oil reservoirs during thermal recovery processes and improve the accuracy and efficiency of temperature transient analysis.