Lake surface water temperature (LSWT) responds rapidly to changes in climatic variables. This response is heterogeneous in space and its spatial distribution is primarily influenced by lake bathymetry and latitude. Such heterogeneity is not captured by one-dimensional water temperature models, which can accurately predict only the average LSWT. We performed a spatially distributed application of the hybrid physically based/data-driven model air2water to predict the LSWT variability in the 5 Laurentian Great Lakes and to deepen our understanding of the role of local depth and latitude in shaping this heterogeneous response. Daily remotely sensed LSWT data were used to calibrate and validate the model during 1995–2018, and additional simulations considering a synthetic warmer climate scenario in which air temperature was increased by 2 °C were run to assess the inter- and intra-lake differences in LSWT warming rates. The model reproduces the observed spatial distribution of LSWT with an average root mean squared error of 1.2 °C and suggests that, under the warmer scenario, the LSWT of the 5 lakes could increase heterogeneously, with the deepest zones showing the maximum warming rates. Summer stratification lengthening is expected to increase with higher local depth; this behaviour attenuates with increasing latitude, whereas the LSWT warming is essentially dependent on the local depth, irrespective of latitude. We highlight the importance of accounting for LSWT spatial heterogeneity to adequately assess the thermal response of the Great Lakes to a warming climate.