Warming trends in the Laurentian Great Lakes and surrounding areas have been observed in recent decades, and concerns continue to rise about the pace and pattern of future climate change over the world's largest freshwater system. To date, most regional climate models used for Great Lakes projections either neglected the lake-atmosphere interactions or are only coupled with a 1-D column lake model to represent the lake hydrodynamics. This study presents a Great Lakes climate change projection that has employed the two-way coupling of a regional climate model with a 3-D lake model (GLARM) to resolve 3-D hydrodynamics essential for large lakes. Using the three carefully selected Coupled Model Intercomparison Project Phase 5 (CMIP5) general circulation models (GCMs), we show that the GLARM ensemble average substantially reduces surface air temperature and precipitation biases of the driving GCM ensemble average in present-day climate simulations. The improvements are not only displayed from an atmospheric perspective but are also evident in the accurate simulations of lake temperature and ice coverage. We further present the GLARM projected climate change for the mid-21st century (2030–2049) and the late 21st century (2080–2099) in the Representative Concentration Pathway (RCP) 4.5 and RCP 8.5 scenarios. Under RCP 8.5, the Great Lakes basin is projected to warm by 1.3–2.1 ^∘C by the mid-21st century and 4.1–5.0 ^∘C by the end of the century relative to the early century (2000–2019). Moderate mitigation (RCP 4.5) reduces the mid-century warming to 0.8–1.8 ^∘C and late-century warming to 1.8–2.7 ^∘C. Annual precipitation in GLARM is projected to increase for the entire basin, varying from 0 % to 13 % during the mid-century and from 9 % to 32 % during the late century in different scenarios and simulations. The most significant increases are projected in spring and fall when current precipitation is highest and a minimal increase in winter when it is lowest. Lake surface temperatures (LSTs) are also projected to increase across the five lakes in all of the simulations, but with strong seasonal and spatial variability. The most significant LST increases occur in Lakes Superior and Ontario. The strongest warming is projected in spring that persists into the summer, resulting from earlier and more intense stratification in the future. In addition, diminishing winter stratification in the future suggests the transition from dimictic lakes to monomictic lakes by the end of the century. In contrast, a relatively smaller increase in LSTs during fall and winter is projected with heat transfer to the deep water due to the strong mixing and energy required for ice melting. Correspondingly, the highest monthly mean ice cover is projected to reduce to 3 %–15 % and 10 %–40 % across the lakes by the end of the century in RCP 8.5 and RCP 4.5, respectively. In the coastal regions, ice duration is projected to decrease by up to 60 d.

The Laurentian Great Lakes are the world's largest freshwater system and regulate the climate of the Great Lakes region, which has been increasingly experiencing climatic, hydrological, and ecological changes. An accurate mechanistic representation of the Great Lakes thermal structure in Regional Climate Models (RCMs) is paramount to studying the climate of this region. Currently, RCMs have primarily represented the Great Lakes through coupled one‐dimensional (1D) column lake models; this approach works well for small inland lakes but is unable to resolve the realistic hydrodynamics of the Great Lakes and leads to inaccurate representations of lake surface temperature (LST) that influence regional climate and weather patterns. This work overcomes this limitation by developing a fully two‐way coupled modeling system using the Weather Research and Forecasting model and a three‐dimensional (3D) hydrodynamic model. The coupled model system resolves the interactive physical processes between the atmosphere, lake, and surrounding watersheds; and validated against a range of observational data. The model is then used to investigate the potential impacts of lake‐atmosphere coupling on the simulated summer LST of Lake Superior. By evaluating the difference between our two‐way coupled modeling system and our observation‐driven modeling system, we find that coupled‐lake atmosphere dynamics can lead to a higher LST during June‐September through higher net surface heat flux entering the lake in June and July and a lower net surface heat flux entering the lake in August and September. The unstratified water in June distributes the entering surface heat flux throughout the water column leading to a minor LST increase, while the stratified waters of July create a conducive thermal structure for the water surface to warm rapidly under the higher incoming surface heat flux. This research provides insight into the coupled modeling system behavior, which is critical for enhancing our predictive understanding of the Great Lakes climate system.

Lake heatwaves (extreme hot water events) can substantially disrupt aquatic ecosystems. Although surface heatwaves are well studied, their vertical structures within lakes remain largely unexplored. Here we analyse the characteristics of subsurface lake heatwaves (extreme hot events occurring below the surface) using a spatiotemporal modelling framework. Our findings reveal that subsurface heatwaves are frequent, often longer lasting but less intense than surface events. Deep-water heatwaves (bottom heatwaves) have increased in frequency (7.2 days decade⁻¹), duration (2.1 days decade⁻¹) and intensity (0.2 °C days decade⁻¹) over the past 40 years. Moreover, vertically compounding heatwaves, where extreme heat occurs simultaneously at the surface and bottom, have risen by 3.3 days decade⁻¹. By the end of the century, changes in heatwave patterns, particularly under high emissions, are projected to intensify. These findings highlight the need for subsurface monitoring to fully understand and predict the ecological impacts of lake heatwaves.

To understand future summer precipitation changes over the Great Lakes Region (GLR), we performed an ensemble of regional climate simulations through the Pseudo‐Global Warming (PGW) approach. We found that different types of convective precipitation respond differently to the PGW signal. Isolated deep convection (IDC), usually concentrated in the southern domain, shows an increase in precipitation to the north of the GLR. Mesoscale convective systems (MCSs), usually concentrated upwind of the GLR, shift to the downwind region with increased precipitation. Thermodynamic variables such as convective available potential energy (CAPE) and convective inhibition energy (CIN) are found to increase across almost the entire studied domain, creating a potential environment more favorable for stronger convection systems and less favorable for weaker ones. Meanwhile, changes in the lifting condensation level (LCL) and level of free convection (LFC) show a strong correlation with variations in convective precipitation, highlighting the significance of these thermodynamic factors in controlling precipitation over the domain. Our results indicate that the decrease in LCL and LCF in areas with increased convective precipitation is mainly due to increased atmospheric moisture. In response to the prescribed warming perturbation, MCSs occur more frequently downwind, while localized IDCs exhibit more intense rain rates, longer durations, and larger rainfall area.

Lake‐effect snow (LES) storms, characterized by heavy convective precipitation downwind of large lakes, pose significant coastal hazards with severe socioeconomic consequences in vulnerable areas. In this study, we investigate how devastating LES storms could evolve in the future by employing a storyline approach, using the LES storm that occurred over Buffalo, New York, in November 2022 as an example. Using a Pseudo‐Global Warming method with a fully three‐dimensional two‐way coupled lake‐land‐atmosphere modeling system at a cloud‐resolving 4 km resolution, we show a 14% increase in storm precipitation under the end‐century warming. This increase in precipitation is accompanied by a transition in the precipitation form from predominantly snowfall to nearly equal parts snowfall and rainfall. Through additional simulations with isolated atmospheric and lake warming, we discerned that the warmer lake contributes to increased storm precipitation through enhanced evaporation while the warmer atmosphere contributes to the increase in the storm's rainfall, at the expense of snowfall. More importantly, this shift from snowfall to rainfall was found to nearly double the area experiencing another winter hazard, Rain‐on‐Snow. Our study provides a plausible future storyline for the Buffalo LES storm, focusing on understanding the intricate interplay between atmospheric and lake warming in shaping the future dynamics of LES storms. It emphasizes the importance of accurately capturing the changing lake‐atmosphere dynamics during LES storms under future warming.

Accurate estimates for the lake surface temperature (LST) of the Great Lakes are critical to understanding the regional climate. Dedicated lake models of various complexity have been used to simulate LST but they suffer from noticeable biases and can be computationally expensive. Additionally, the available historical LST datasets are limited by either short temporal coverage (<30 years) or lower spatial resolution (0.25° × 0.25°). Therefore, in this study, we employed a deep learning model based on Long Short-Term Memory (LSTM) neural networks to produce a daily LST dataset for the Great Lakes that spans an unparalleled 42 years (1979–2020) at a spatial resolution of ~1 km. In our dataset, the Great Lakes are represented by ~33,000 unstructured grid points and the LSTM training incorporated the information from each grid point. The LSTM was trained with seven meteorological variables from reanalysis data as feature variables and the LST from a historical satellite-derived dataset as the target variable. The LSTM was able to capture the spatial heterogeneity of LST in the Great Lakes well and exhibited high correlation (≥0.92) and low bias (limited to ±1.5 °C) for the temporal evolution of LST during the training (1995–2020) and testing (1979–1994) periods.