The effect of natural convection on the temperature distribution of skin was numerically investigated by solving the conjugate heat transfer of the skin under contrast therapy, where hot and cold stimulations are periodically induced on the skin. A finite volume method based on the SIMPLE algorithm was adopted to solve the axisymmetric incompressible Navier-Stokes equations coupled with an energy equation. Those equations were strongly coupled with the Pennes bio-heat equation of skin for the analysis of the conjugate heat transfer. The amplitudes of the sinusoidal temperature profiles of the skin at selected depths were found to become smaller when natural convection was considered, and the temperature evolution obtained by analyzing the conjugate heat transfer differed from that obtained with constant heat transfer coefficients or with a conduction effect only, especially when blood perfusion was included. The spatially-averaged heat transfer coefficient predicted from the present simulation of conjugate heat transfer of contrast therapy was also confirmed to be close to the typical values employed in previous numerical studies. More importantly, the present conjugate heat transfer simulation for contrast therapy revealed that the heat transfer coefficient on skin varies both spatially and temporally. Therefore, the conjugate heat transfer analysis, in which the heat transfer coefficient is calculated as a function of time and space, must be employed in order to accurately predict the temperature evolution inside the skin when subjected to contrast therapy. (C) 2015 Elsevier Ltd. All rights reserved.