Abstract While changing lake surface conditions have received significant research scrutiny, changes in subsurface conditions, including stratification and heat content, remain largely unexplored. In this work, we highlight changes in thermal structure, stratification dynamics, and ice characteristics in the Laurentian Great Lakes (Lakes Superior, Huron, Michigan, Erie, and Ontario) as simulated between 1979 and 2021. Three-dimensional lake hydrodynamics and ice cover were modeled using the Finite Volume Community Ocean Model (FVCOM) coupled with the Los Alamos sea ice model (CICE). Analysis revealed significant increases in surface (0.4°–0.6°C decade−1) and subsurface (0.1°–0.4°C decade−1) temperatures as well as dramatic losses in ice cover (1%–8% decade−1) and ice volume (0–3 km3 decade−1) over the last 40 years. Estimated surface heating rates were strongest during the summer and fall, while subsurface warming was most rapid during the nearly isothermal winter and spring. Intensified (decreased) summer (winter) stratification led to shifts in lake turnover dynamics, with delayed fall turnover dates (2–6 days decade−1) and earlier spring overturn dates (2–9 days decade−1). Modeled surface temperatures (LST), bottom temperatures (LBT), and annual averaged ice cover (AAIC) were used to estimate low-frequency climate signals, which were highly correlated with the Atlantic multidecadal oscillation. Warming trends fit to residual climate signals (LST: 0.1°C decade−1; LBT: 0.03°C decade−1; AAIC: −1% decade−1), calculated by removing low-frequency variability from the raw climate signal, were lower than those fit to associated low-frequency components, suggesting that recent climate change in the Great Lakes may be strongly influenced by natural multidecadal climate variability.