Temperature distribution in subduction zones is a key factor in understanding magma generation, earthquake occurrence, metamorphic reactions, and geochemical element recycling. In this thesis, I have developed a new numerical method to model subduction zone thermal structures with real slab geometries. The new method couples the finite-difference method with the finite-element method to solve the conduction-convection heat transfer system involved in a subduction zone. The mantle and wedge convection is simulated with the finite-element method with the incoming slab convergence rate imposed as a boundary condition at the slab surface. The heat transfer problem is solved with the finite-difference method. It uses a staggered-grid discretization approach so that the effect of the mantle and wedge convection on the thermal model can be accurately taken into account. With material averaging, the simulation method can easily handle a curved slab with a simple, structured grid. The thermal structures that I have calculated for the east Aleutian subduction zone show that a two-segment slab model can overestimate the slab surface temperature at 100 km by up to 90°C for a 75-km thick overriding lithosphere. The thermal structures of ten subduction zones around the Pacific Rim have been generated including the east Aleutians, the Cascades, the Central and South America, the northeast Japan, and Mariana subduction zones. I have shown that the predicted slab surface temperature at 100 km is correlated with B/Zr ratio. As the slab surface temperature increases, the B/Zr ratio decreases systematically for two different levels of B-enrichment (5 ppm and 10 ppm). The results show that the slab tip temperatures have little correlation with the slab lengths and depths.