When converting temperature-time curves obtained from geochronology into the denudation history of an area, variations in the isotherm geometry should not be neglected. The geothermal gradient changes with depth due to heat production and evolves with time due to heat advection, if the deformation rate is high. Furthermore, lateral variations arise due to topographic effects. Ignoring these aspects can result in significant errors when estimating denudation rates. We present a numerical model for the thermal response to thrust faulting, which takes these features into account. This kinematic two-dimensional model is fully time-dependent, and includes the effects of alternating fault activation in the upper crust. Furthermore, any denudation history can be imposed, implying that erosion and rock uplift can be studied independently to each other. The model is used to investigate the difference in thermal response between scenarios with simultaneous compressional faulting and erosion, and scenarios with a time lag between rock uplift and denudation. Hereby, we aim to contribute to the analysis of the mutual interaction between mountain growth and surface processes. We show that rock uplift occurring before the onset of erosion might cause 10% to more than 50% of the total amount of cooling. We applied the model to study the Cenozoic development of the Sierra de Guadarrama in the Spanish Central System, aiming to find the source of a cooling event in the Pliocene in this region. As shown by our modeling, this temperature drop cannot be caused by erosion of a previously uplifted mountain chain: the only scenarios giving results compatible with the observations are those incorporating active compressional deformation during the Pliocene, which is consistent with the ongoing NW-SE oriented convergence between Africa and Iberia. © 2004 Elsevier B.V. All rights reserved.