In the present study a parallel three dimensional Volume of Fluid (VOF) method is developed to simulate Marangoni convection in immiscible fluids with variable surface tension. Conservation equations are solved based on cell-averaged one-field volume tracking scheme. Evaluating the convective term in the energy equation along the boundary between the fluids highly depends on the position and orientation of the interface; hence, using average cell values simply ignores the interface shape and leads to computational uncertainty. As a remedy to this issue, the original idea behind the volume tracking method is used not only to advect mass and momentum but also energy across cells. To verify the proposed algorithm, results are compared against theoretically predicted thermocapillary migration velocity of a droplet at the limit of zero Marangoni number. However, at relatively high Marangoni numbers, thermal boundary layers are very thin and challenging to resolve. To demonstrate the capabilities of the heat transfer module, simulation of a Fluorinert droplet moving in silicon oil under applied temperature gradient is compared against the available experimental results and the migration velocity of the droplet in microgravity is reported. (C) 2014 Elsevier Ltd. All rights reserved.