Rapidly expanding lake heatwaves under climate change

Lake heatwaves—prolonged periods of hot surface water temperature in lakes—have recently been shown to increase in intensity and duration, with numerous potential implications for aquatic ecosystems. However, an important physical attribute of lake heatwaves that has not yet been investigated is their spatial extent, and how it varies within a warming world. Here, we show that the spatial extent of lake heatwaves, defined as contiguous regions within a lake that simultaneously experience extreme warm conditions, is increasing in the largest group of freshwater lakes on Earth, The Laurentian Great Lakes. We show that the maximum spatial extent of lake heatwaves is sensitive to inter-annual variations in winter ice cover and the timing of stratification onset in spring. Notably, we find that a lengthening of the warm summer season and, in turn, an overall increase in surface water temperature, stimulates the development of larger lake heatwaves. On average, our results suggest that the mean spatial extent of lake heatwaves has increased two-fold since 1995. We anticipate this rapid expansion of lake heatwaves to have widespread implications for heat-related impacts on aquatic species.


Introduction
Much of the focus of global warming impacts on lakes have concentrated on seasonal or annual changes in mean lake surface water temperature, with previous studies suggesting a long-term warming of lakes worldwide (Schneider and Hook 2010, Zhang et al 2014, O'Reilly et al 2015, Woolway and Maberly 2020, Dokulil et al 2021. In comparison, little is known about lake heatwaves-periods of warm lake surface water temperature-and how they respond to climatic variations. One notable exception is the study of Woolway et al (2021a) which, using a onedimensional lake model, suggested that lake heatwaves become more intense and longer lasting within a warming world. The rise of unprecedented temperatures, as occur during lake heatwaves, represents a potential trigger for catastrophic mass mortality events, such as fish die-offs, as lake temperatures exceed organismal physiological thresholds (Till et al 2019). Over time, an increase in the exposure of lake ecosystems to extreme temperatures may lead to an irreversible loss of species, as has already been observed in the oceans due to the increased frequency of marine heatwaves , Straub et al 2019. Lake heatwaves could also promote the growth of harmful cyanobacteria, which are adapted to endure warmer and irregular temperatures (Jöhnk et al 2008, Duan et al 2009, Zhang et al 2016, Rasconi et al 2017, Baker and Geider 2021, thus having negative impacts on some of the key benefits that lakes provide to society, such as the provision of safe water for drinking and irrigation. Given the vulnerability of aquatic ecosystems to thermal extremes, an understanding of lake heatwaves has taken on new and critical importance.
An important physical attribute of lake heatwaves that has not previously been investigated is their spatial extent, specifically, the spatial extent of contiguous regions within a lake that simultaneously experience heatwave conditions. The physical size of such regions has important implications for heat-related impacts on aquatic species, including their exposure to extreme warm water temperature. Here, we address this critical knowledge gap by investigating the spatial extent of heatwaves in the largest group of freshwater lakes on Earth, The Laurentian Great Lakes. In this study we calculate the spatial extent of lake heatwaves observed each year since 1995 and investigate how they respond to climatic variations. We also relate these observed changes to the inter-annual variability in winter ice cover and the onset of thermal stratification (i.e. the time of year in which the positively stratified season starts), two key metrics known to influence the thermal response of The Laurentian Great Lakes to climate change (Austin and Colman 2007, Anderson et al 2021.

Study system
The Laurentian Great Lakes (situated in the upper mid-east region of North America; figure S1 (available online at stacks.iop.org/ERL/16/094013/ mmedia)), hereafter referred to as the Great Lakes, are the largest freshwater system on Earth with respect to surface area and contain approximately 21% of the world's unfrozen surface freshwater. Lake Michigan-Huron, when considered as a single system, and Lake Superior make up the two largest lakes by surface area, but as a whole the Great Lakes contain over 7000 km of coastline that stretches across the border between Canada and the United States. Lake Superior is the deepest of the lakes, with regards to average depth (=149 m), followed by Ontario (=86 m), Michigan (=85 m), Huron (=59 m), and Erie (=19 m). The Great Lakes are generally considered dimictic, mixing in the fall and spring, with a high degree of interannual and inter-lake variability in winter ice cover. Each lake has experienced a decreasing longterm trend in ice cover since the 1970s. Between two Canadian provinces and eight US states, the Great Lakes system provides drinking water to over 30 million people and supports a 7-billion-dollar fishery. Changes to the physical characteristics of the lakes, in particular temperature and ice, can have significant impacts to water quality, the ecosystem, and economy.

Surface water temperature observations
Spatially resolved daily lake surface water temperatures were acquired from the Great Lakes surface environmental analysis (GLSEA) for the period 1995-2020 (Schwab et al 1999). This data set is commonly used for investigating long-term surface temperature trends in the Great Lakes (Mason et al 2016, Anderson et al 2021. The GLSEA is a ∼1.8 km resolution satellite-derived surface temperature product that provides uninterrupted time series of daily lake surface temperatures. The data is freely available at https://coastwatch.glerl.noaa.gov/glsea/. The satellite data set is based on the NOAA advanced very highresolution radiometer, a sensor aboard NOAA's Polar Orbiting Environmental Satellites. As with other satellite data products, the surface water temperatures are only derived from satellite remote sensing during cloud/rain free periods. Notably, the GLSEA product uses cloud/rain free pixels to derive surface water temperature, and then uses a spatiotemporal interpolation process to provide complete daily fields of surface temperature for each lake. These surface temperature fields are continually validated against in-situ (buoy) observed water temperatures. Most notably, NOAA and Environment and Climate Change Canada (ECCC) operate several buoys in the Great Lakes (www.ndbc.noaa.gov/) that provide realtime surface water temperature. The GLSEA product was developed using these observed water temperatures and now as an operational NOAA product it is continually (daily-tracked) verified against these surface measurements.

Lake heatwaves
Using the GLSEA temperature data, we identified lake heatwaves for each 1.8 km grid following the methods described by Hobday et al (2016), and as used by Woolway et al (2021a). Specifically, using the R package 'heatwaveR' (Schlegel and Smit 2018), lake heatwaves were defined as when daily mean lake surface temperatures were above a local and seasonally varying 90th percentile threshold for at least five consecutive days (figures S2 and S3). Two lake heatwave events with a break of less than three days were considered as a single event. Lake surface temperature anomalies were calculated for each calendar day using the daily temperatures within an 11 d window centered on the date across all years within a defined period and smoothed by applying a 31 d moving average (Hobday et al 2016). For each day, a binary map that identified which of the 1.8 km grids across the Great Lakes experienced a heatwave was generated (figure S4). The 3D binary maps (longitude × latitude × time) were then used to identify spatiotemporally connected heatwaves (i.e. grids that are connected in 3D space in terms of heatwave occurrence) following a connected component analysis. Specifically, a connected component analysis was used to identify a connected group of heatwave pixels which, in turn, were part of the same heatwave event i.e. a contiguous region simultaneously experiencing a heatwave (figure S5). For each lake, we then identify the largest heatwave that occurred each year, and calculated how the spatial extent of these heatwaves varied from 1995 to 2020. Note, that in this study we focus only on the largest spatiotemporally connected heatwaves, as these are the events that are likely to have the largest impact on the aquatic ecosystem. Smaller and isolated lake heatwave events in each lake, such as those occurring in an individual grid or a small number of connected grids, were not considered in this study. Also, we note that as the satellite data are not available during periods of cloud/rain-at which time a  (a)) spatial extent of the largest heatwaves; and (d) the spatial coverage of the largest heatwave in each lake at the time of their maximum extent, with the colors representing the heatwave intensity (i.e. the temperature anomaly relative to the long-term average at that time of year). The gray regions in panel (d) illustrate that the heatwave did not reach that region of the lake. The spatial extent of lake heatwaves is shown as percentages, calculated relative to the total surface area of each lake.
smoothing algorithm is applied-our study does not strictly consider changes in lake thermal extremes during storm events.

Drivers of lake heatwave extent
The spatial extent of the largest heatwaves in each lake were compared with inter-annual variations in the annual maximum ice cover percentage, the date of stratification onset, and the over-lake air temperatures. Annual maximum ice cover for the Great Lakes (1995-2020) were acquired from the Great Lakes Environmental Research Laboratory (www.glerl.noaa.gov/data/ice/) and represent the maximum surface area covered by lake ice each year, and is a useful proxy for summarizing winter conditions in the Great Lakes. The average date of stratification onset per lake were calculated from the lake-average GLSEA data, as the first day of year in which the surface water temperatures exceeded the temperature of maximum density of freshwater (3.98 • C), as defined in previous studies (Austin and Colman 2008, Woolway and Merchant 2017, Fichot et al 2019, Woolway et al 2021b. Air temperature observations for each of the Great Lakes were acquired from the Great Lakes Coastal Forecast System (GLCFS) (Schwab and Bedford 1994). This data contains observations from several over-water and coastal meteorological stations that include NOAA Coastal-Marine Automated Network (C-MAN) sites, National Weather Service buoys, the Automated Surface Observing System, and United States Coast Guard Stations (NOAA 2020). In this study, we use sites in close proximity to each of the Great Lakes (figure S1), in which the GLCFS applies atmospheric adjustments to a common measurement height and employs a natural neighbor interpolation (Sambridge et al 1995) to create hourly overwater meteorological conditions on a 2 km grid, which are then spatially averaged to provide lake-averaged air temperatures. Throughout the analysis, anomalies are calculated for each variable relative to the average conditions for the day/month of year evaluated over the 1995-2020 base period.

Results
To illustrate the approach followed in this study for calculating the long-term change in the spatial extent of lake heatwaves we show, in figure 1, the computed heatwave spatial extent in each studied site during the first year of observations. The spatial extent of the largest lake heatwaves in 1995, calculated as a percentage of the total lake surface area, varied temporally from its inception to its end ( figure 1(a)). Likewise, the computed mean and maximum spatial extent of the largest heatwaves observed in each lake varied considerably. The largest lake heatwave in 1995 was observed in Lake Erie, which occupied a region equivalent to 94% of the total lake surface area at its peak ( figure 1(c)). On average, from its start to end, this large heatwave occupied a region equivalent to 37% of the total lake surface area ( figure 1(b)). By repeating the approach described above for each year with observational data we explore, for all lakes, the inter-annual variability in the mean and maximum spatial extent of lake heatwaves. As Colors illustrate the heatwave intensity (i.e. temperature anomaly). Gray regions illustrate that the lake heatwave did not reach that region of the lake. Table 1. The largest heatwaves observed in the Great Lakes since 1995. Shown are the maximum percentage coverage of the largest heatwaves observed in each of the Great Lakes. The spatial extent of lake heatwaves is shown as percentages, calculated relative to the total surface area of each lake. The maximum heatwave percentages are shown for the five years with the largest events and ranked according to their spatial extent. All years are shown in table S1.

Superior
Michigan  ; table S1). In other years, lake heatwaves can occupy a vast proportion of the total lake surface area. For example, in 2012, approximately 99% of the surface area of Lakes Superior, Michigan and Ontario was occupied by a lake heatwave at its peak extent (figure 2; table 1). On average, across the Great Lakes, ∼98% of the total surface area experienced a heatwave in 2012 (table 1). According to their maximum spatial extent, five of the largest lake heatwaves averaged across the Great Lakes have occurred since 2010 (table 1). However, some exceptions can be noted in individual lakes, such as the anomalously large (covering 92% of the lake surface area) heatwave observed in Lake Superior in 1998 (figure S6; table 1).
To explore the mechanistic drivers of lake heatwaves, we now compare the inter-annual variability of their maximum spatial extent, averaged across the Great Lakes, to the inter-annual variability in the percentage of maximum winter ice cover and the subsequent date of stratification onset (figure 3). Here, we observe a statistically significant (r 2 = 0.3; p = 0.003) negative relationship between the interannual variability in the maximum heatwave extent and the percentage of winter ice cover ( figure 3(a)). Specifically, during years when ice cover is low, the maximum spatial extent of lake heatwaves is high. This association is driven by the fact that less winter ice cover, in addition to warmer surface air temperatures in spring, contributes to an earlier onset of thermal stratification ( figure 3(b)). An earlier stratification onset, as well as higher surface air temperatures during the open-water season, likewise promotes warmer surface water temperature (figure 3(c)). Higher surface water temperatures, in turn, stimulate the development of large lake heatwaves (figure 3(d)). We also find that the deepest lakes experience the greatest association between stratification onset and the maximum spatial extent of lake heatwaves, evaluated via a linear regression model between the interannual variability in the date of stratification onset and the heatwave maximum spatial extent. Notably, the proportionate amount of variation in the maximum spatial extent explained by the onset of stratification is lowest in the shallowest lake (Erie) and highest in the deepest lake (Superior) (figure S9).
The average spatial extent of lake heatwaves across the Great Lakes demonstrates a noticeable increase since 1995, albeit with considerable inter-annual variability ( figure 4(a)). Most notably, our observations suggest that the spatial extent of lake heatwaves in the Great Lakes has increased two-fold between the first and last decade (1995-2004 and 2011-2020, respectively) of the study period, from 10.8% to 26.0% on average. However, we also highlight that the average spatial extent is much greater in some specific years. For example, in 2016 the average spatial extent of lake heatwaves across the Great Lakes was 62%, meaning that almost two thirds of the Great Lakes' surface area were occupied by a heatwave. Here, we also calculate that the rate of change in the average spatial extent of lake heatwaves across the Great Lakes is 7.3% decade −1 . The rate of change is greatest in the shallowest of the Great Lakes (Erie) at 11.4% decade −1 , and is least in the deepest lake (Superior) at 3.5% decade −1 . While the number of studied lakes included in this investigation does not allow for a comprehensive analysis of this depth dependence, we do find a clear relationship between the rate of change in the average spatial extent of lake heatwaves and depth of the Great Lakes (figure S10). Intuitively, the average spatial extent of the observed heatwaves across the Great Lakes is also strongly related to their average duration ( figure 4(b)). Thus, as lake heatwaves become increasingly long-lived within a warming world, their average spatial extent will likely increase further. Inter-annual variability in the mean spatial extent of lake heatwaves. Shown for each lake, as well as the average across the Great Lakes, is (a) the inter-annual variability in the mean spatial extent of lake heatwaves, calculated as the average heatwave coverage over the duration of the heatwave event. We also show (b) the relationship between the mean spatial extent of lake heatwaves (i.e. across the Great Lakes) and the average duration of these lake heatwaves. The spatial extent of lake heatwaves is shown as percentages, calculated relative to the total surface area of each lake.

Discussion
For the first time, we investigated the spatial extent of lake heatwaves and how they respond to climate change. Our analysis suggests that the mean spatial extent of heatwaves in the largest group of freshwater lakes on Earth has increased since 1995, by up to two-fold in the most recent decade. This increase in average spatial extent is driven primarily by an increase in the duration of lake heatwaves which, ultimately, allows them to grow to unprecedented levels in recent years. During some extreme warm years, lake heatwaves can occupy almost the entire surface area of the Great Lakes for a sustained period but, in contrast, can cover a relatively small area during other, notably cooler, years. The inter-annual variability in the maximum spatial extent of lake heatwaves is driven primarily by the inter-annual variability in the date of stratification onset as well as the magnitude of air temperature anomalies during the open-water season. An earlier onset of stratification, which can occur following mild winters and/or anomalously warm spring/summer air temperatures (Woolway et al 2021b), can result in a lengthening of the summer stratified season and, in turn, warmer surface water temperatures (Austin and Colman 2007, Piccolroaz et al 2015, Zhong et al 2016, Woolway and Merchant 2017. Most notably, the rate of surface warming in seasonally stratifying lakes is substantially greater after the onset of thermal stratification when the volume of water that participates directly in air-water surface heat exchange is smaller. Thus, an earlier stratification onset in lakes can result in an increase in the time during which the shallow upper mixed layer is present and, in turn, result in warmer lake surface temperatures which, subsequently, promotes the development of larger lake heatwaves. This is clearly evident in Lake Superior in 1998, where warmer than average surface water temperatures driven by the aforementioned mechanisms (Van Cleave et al 2014, Piccolroaz et al 2015, Zhong et al 2016, contributed to the development of an anomalously large lake heatwave. Our analysis also identified a depth effect on the association between stratification onset and the maximum spatial extent of lake heatwaves. The maximum spatial extent of lake heatwaves is more strongly correlated with the date of stratification onset in the deepest lakes, which is due to temperature anomalies associated with stratification onset persisting for longer in deeper lakes where thermal inertia is large (Toffolon et al 2014, Woolway andMerchant 2017).
We expect other seasonally stratifying lakes globally to experience a similar increase in the spatial extent of heatwaves in recent decades, if the mechanistic drivers described above, notably an earlier onset of stratification (Woolway et al 2021b), or an increase in surface water temperature due to other atmospheric of within-lake drivers , are changing in line with those observed in our studied sites. Ultimately, an increase in lake surface water temperature will result in the occurrence of longer lasting lake heatwaves, which will result in an increase in their spatial extent. In lakes that mix frequently, at daily to weekly intervals (i.e. polymictic lakes), which according to Hutchinson's (1957) mixing classification are among the most abundant at mid-latitudes, it is unlikely that the spatial extent of heatwaves will increase within a warming world, as vertical mixing will likely destroy any spatiotemporally connected heatwaves that develop. In lakes that are semi-permanently/permanently stratified (oligomictic/meromictic), such as the deep lakes south of the European Alps or those located at low latitudes (e.g. Lake Tanganyika), it is possible that lake heatwaves could respond sensitively to a warming world, with lake heatwaves extending across multiple years, as has previously been observed in the marine environment (Bond et al 2015, Di Lorenzo and Mantua 2016, Oliver et al 2021. Moreover, as lakes transition to a more stable mixing regime in the future (e.g. from polymictic to monomictic or monomictic to oligomictic; Woolway and Merchant 2019), more lakes will likely experience an expansion in the spatial extent of heatwaves. This is an area for further study.
We anticipate the spatial extent of lake heatwaves to increase at a faster rate than observed here in the Great Lakes. This is particularly true in lakes that will lose their ice cover this century (Sharma et al 2019), as lake heatwaves could then extend throughout the winter season and continue to expand. Based on our results, one might also expect that, in the near future, lake heatwaves may occupy the entire surface area of lakes on an annual basis, potentially resulting in the emergence of a permanent heatwave state (Woolway et al 2021a). The continued spatial expansion of lake heatwaves will likely expose aquatic organisms to thermal extremes that have not been experienced before, and there will likely be a substantial departure from the 'normal' thermal conditions that have shaped the aquatic environment in the past. Aquatic organisms that live in regions close to their maximum thermal limit are likely to be most susceptible to the rapid expansion of lake heatwaves, as temperature anomalies there are more likely to exceed organismal physiological thresholds and hence increase mortality rates (Till et al 2019). Specifically, with the spatial expansion of lake heatwaves, one could expect dramatic losses of marginal populations at thermal range edges in the Great Lakes, as has already been documented in the marine environment (Hobday et al 2016). A potential consequence of this rapid expansion is a reduction of habitat for cold-waters species, which may also result in a compositional shift to warmer-water species that are potentially non-native (Cline et al 2013, Elliott et al 2015. However, gaps in our current understanding of how lake heatwaves will affect aquatic organisms and communities limit our ability to forecast the precise effects of this anticipated spatial expansion, notably due to the complex interactions that occur within a lake in response to climate change. In addition, during summer stratification, cooler water at depth could provide a potential thermal refuge (Kraemer et al 2021), thus allowing aquatic organisms to escape the extreme conditions associated with lake heatwaves. However, we must also consider that environmental conditions in the deeper waters during a lake heatwave event may not always be suitable with, for example, many lakes worldwide experiencing oxygen depletion at depth during the summer season (Jane et al 2021). While the exact negative impacts of lake heatwaves under climate change is unknown, it is likely that when combined with the host of other human effects, such as eutrophication, aquatic ecosystems will be adversely affected by the spatial expansion of lake heatwaves.

Conclusions
We analyzed spatially resolved surface water temperatures from the Great Lakes to investigate how the spatial extent of lake heatwaves-defined as contiguous regions within a lake that simultaneously experience extreme warm conditions-vary under climate change. Our analysis revealed that the maximum spatial extent of lake heatwaves is sensitive to interannual variations in winter ice cover and the timing of stratification onset. Notably, our findings suggest that a lengthening of the warm summer season and, in turn, an overall increase in surface water temperature, results in the development of larger spatiotemporally connected heatwaves. Furthermore, our analysis suggests that the mean spatial extent of lake heatwaves has increased two-fold, from 10.8% to 26.0% on average, between the first and last decade of the study period. The rate at which the mean spatial extent of lake heatwaves is increasing varies across the Great Lakes from 3.5% decade −1 in Lake Superior to 11.4% decade −1 in Lake Erie. On average, the spatial extent of lake heatwaves is increasing at a rate of 7.3% decade −1 . Given the sensitivity of lake ecosystems to hot temperature extremes, we anticipate this observed expansion of lake heatwaves to have widespread implications for heat-related impacts on aquatic species.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: http:// doi.org/10.5281/zenodo.5115909.