The presence of anthropogenic infrastructure within the ground under cities can act to increase the temperature of the subsurface in a phenomenon known as the subsurface urban heat island effect. Large city-scale numerical subsurface temperature models accounting for both the heat fluxes stemming from human-made structures as well as hydrogeological variation are often computationally prohibitive due to the vastly different scales of modelling required. We present a novel and scalable methodology for the determination of the subsurface thermal climate beneath urban centres at a holistic scale, grouping volumes of ground into archetypes according to common thermal behaviours, identified using supervised machine learning algorithms. We apply this methodology to two boroughs in the centre of London, UK, with the results showing good agreement with a higher-resolution model for one of the boroughs. We further demonstrate the scalability of the approach by modelling the boroughs sequentially, where the second borough was modelled at little to no additional computation cost at a loss of accuracy in the local depth -averaged temperature between 0.10 and 0.23 degrees Celsius. This methodology is an important step towards determining large-scale subsurface climate beneath cities, providing insights for various applications, such as an improved understanding of geothermal energy potential.