We study crustal thermal evolution by examining heat flow patterns along a convergent boundary from a young subduction zone to a more structurally mature collision zone. More than 8000 km of seismic profiles covering an offshore region of 45,000 km2 in southern Taiwan show widespread bottom-simulating reflectors (BSRs). We derived 1107 BSR-based heat flows before combining 42 additional, published, offshore thermal probe data and 86 on-land heat flow data to document the shallow forearc thermal structures from the subduction zone to the collision zone. In the subduction zone, the geothermal gradient ranges mostly within 30–80 °C/km, and decreases toward the arc due to slab cooling, intensive dewatering at the toe, sediment blanketing, topographic effects, and other processes. The geothermal gradient ranges mostly from 30 to 90 °C/km in the collision zone, and increases, instead of decreases, toward the arc, possibly caused by exhumation, erosion, topographically induced ground-water circulation, and some upper mantle processes related to collision. Heat flow in the collision zone ranges from 80 to 250 mW/m2. The high heat flow in the collision zone correlates with a shallower seismicity zone and high seismic attenuation, while the lower heat flow in the subduction zone might allow the earthquakes to rupture to greater depth. The heat flow increases along the topographic high from subduction to collision zone due to increasing geothermal gradients and higher thermal conductivities of the exhumed basement rocks. This heat flow variation may generate an artificial exhumation rate pattern, if a conventional 30 °C/km geothermal gradient was used in fission-track studies.