Specifying exact geometry of vessel network and its effect on temperature distribution in living tissues is one of the most complicated problems of bioheat field. Presence of blood and its thermal role in the living tissues such as blood perfusion, convection, perfusion of arterial-venous blood through capillaries and physiological processes, including metabolic heat generation and also heat conduction and interaction with environment, make this process complicated. Therefore in this project the effects of blood vessels on temperature distribution in a skin tissue subjected to various thermal therapy conditions are investigated. A precise and realistic model for skin and blood vessels is used to analyze the transient heat transfer through skin. Counter-current vessels are embedded in a three-dimensional Skin structure including three layers. Branching angles of vessels are calculated using the physiological principle of minimum work. Length and diameter ratios are specified using length doubling rule and Cube law, respectively. Temperature distributions in three-dimensional space are evaluated by solving the continuity, momentum and energy equations for blood flow and modified Pennes bioheat (TWMBT) equation for tissue. TWMBT model assumes heat is propagated by waves with finite speed and therefore there is a relaxation time for temperature elevation in the skin tissue. Since blood reaches to equilibrium with tissue within capillary bed, convection of blood can be neglected. Hence, thermal effects of the capillary bed are replaced with an energy source. Blood is considered non-Newtonian fluid using the power law model. Results for temperature distribution show low temperature contours along the artery meaning that arteries act as heat sinks; while high temperature contours along the veins indicate that veins carry heat out of the upper part of the tissue. Due to low thermal conductivity of biological tissue, lower skin levels are warmed sooner by convection than by conduction. Also, cooling effect of blood in the capillary bed is calculated precisely. It is observed that as the tissue temperature increases, perfusion rate in lowering the tissue temperature enhances. Totally, it is observed that the role of blood for cooling the tissue at high temperature is significant. The effects of different parameters such as vascular network, boundary conditions, relaxation time, thermal properties of skin, metabolism, convection coefficient, ambient temperature and pulse heat flux on temperature distribution are investigated. Tremendous effect of boundary condition type at the lower boundary is noted. It seems that neither insulation nor constant temperature at this boundary can completely describe the real physical phenomena. It is expected that real temperature at the lower levels be somewhat between two predicted values. The effect of temperature on the thermal properties of skin tissue is considered. It is shown that considering temperature dependent values for thermal conductivity is important in the temperature distribution estimation of skin tissue; however the effect of temperature dependent values for heat capacity is negligible. It is seen that considering modified Pennes equation in processes with high heat flux during low time (such as some thermal treatments) is very significant and leads to results completely different from the classic Pennes equation. Keywords: Triple layered skin, counter-current vessels, thermal therapy, modified Pennes equation, relaxation time