We present a new method to incorporate paleo-surface temperature effects in steady state 3D conductive temperature models. The workflow approximates the transient effects and incorporates these into steady state models, using appropriate source and sink terms for radiogenic heat production. This allows for rapid models, which can be easily used in ensemble approaches for data assimilation of high-resolution temperature models for geothermal resource assessment. The workflow is demonstrated for the Netherlands which is a sedimentary basin with a wealth of deep (temperature) data from groundwater and oil and gas wells and past studies on the 3D temperature distribution. 3D subsurface temperature models of the Netherlands ranging up to 5 km depth systematically overpredict temperatures at shallow (<1500 m) depth. Analysis of both shallow (<600 m, >200,000 measurements) and deep temperature measurements (1–6 km, >1500 measurements), clearly demonstrates a shallow thermal gradient in over 200 locations of ~20 °C km−1 for the top 400 m underlain with a deep geothermal gradient of ~31 °C km−1 for the 2–4 km interval. Improvements in 3D subsurface modelling regarding the shallow part are accomplished by adding a paleo-surface temperature correction related to glaciation effects of the Weichselian glacial period. This paleo-surface temperature correction proved to be the missing link between two distinctive geothermal gradients observed and is consistent with earlier findings for limited datasets. The consistent overprediction of modelled temperatures in 74% of locations for the top 2 km which are regularly distributed over the Netherlands demonstrates that the influence of paleo-surface temperatures is rather uniform over large areas and not significantly overprinted by other effects such as groundwater flow. The updated model, marked by up to 10 degrees cooling compared to models ignoring the paleo-surface temperature effects, has major implications for assessing geothermal resource potential up to 2 km ...