Various datasets were used in this study for data analysis, model input, and model validation. The data for SWAN model inputs included hourly surface wind forcing and daily ice cover. The hourly surface wind forcing was obtained from the Climate Forecast System Reanalysis (CFSR)⁴ dataset with a spatial grid resolution of ∼0.3 degree for 1979–2010 and the upgraded Climate Forecast System Version 2 (CFSv2)⁵ dataset with a higher spatial resolution of ∼0.2 degree, which is archived as an extension of CFSR and available from 2011 to present at the National Centers for Environmental Prediction (NCEP) (Saha et al., 2010, 2011). CFSR is a global, high-resolution, coupled atmosphere-ocean-land surface-sea ice system, with the assimilation of satellite radiances and all available conventional and satellite observations. Analysis of the CFSR output indicates the product is far superior in most respects to the reanalysis of the mid-1990s (Saha et al., 2010, 2011), and has shown good performance over the Great Lakes regions (Jensen et al., 2012; Xue et al., 2015). The daily ice cover on Lake Michigan was obtained from the Great Lakes Ice Cover Database,⁶ which includes the Great Lakes Ice Atlas⁷ for the period 1973–2002, with addendums for 2003 through present.

Both the daily ice cover data described above and the monthly lake-wide average water levels were analyzed to understand their impacts on the wave climate. Monthly lake-wide average water levels were obtained from the Great Lakes dashboard of NOAA Great Lakes Environmental Research Laboratory.⁸ These lake-wide averages were derived based on a selected set of U.S. and Canadian station data as determined by the Coordinating Committee on Great Lakes Basic Hydraulic and Hydrologic Data, which represents the best estimates of the lake-wide mean water levels.

Lake-wide buoy observations were used for the SWAN model validation. Surface conditions across Lake Michigan were well documented with two offshore buoys (MetOcean buoys operated by the NDBC of NOAA) and sixteen coastal buoys (operated through the Great Lakes Observing System (GLOS) with funding from the Integrated Ocean Observing System (IOOS) of NOAA). The two offshore MetOcean buoys are NDBC S45002 and S45007. S45002 is a 3-meter discus buoy with Self-Contained Ocean Observing Payload (SCOOP), which provides wind speed, wind direction, wind gusts, wave height, dominant wave period, average wave period, mean wave direction, and a variety of other air/sea measurands hourly. NDBC S45007 is a 2.3-m foam discus buoy also with SCOOP and similarly reports hourly. The sixteen coastal buoys operated by a variety of universities and other GLOS partners are primarily TIDAS 900 buoys.⁹ These buoys consist of a 1.1 m-diameter conical hull with a six-degree-of-freedom wave sensor providing at a minimum, the same suite of sensors as the NDBC buoys, however, at a much higher rate of reporting to capture the rapid change characteristics of coastal regions. Equipped with the SeaView Systems 603 wave sensor,¹⁰ TIDAS buoy motions were recorded for 17.06 min at a sample rate of 2 Hz for a time series sample of 2048 observations. Buoy wave parameters were reported at 20-min intervals and available through GLOS and NDBC.

Prevalent numerical models for wave hindcast are the third-generation spectral wave models including Wave Modeling (WAM, WAMDI Group, 1988), WAVE-WATCH III (Tolman, 1991), TOMAWAC (Benoit et al., 1997) and Simulating Wave Nearshore (SWAN, Booij et al., 1999). In this study, the SWAN model was used to investigate the spatiotemporal changes of wave climate. SWAN is a third-generation spectral wave model developed at Delft University of Technology that computes random, short-crested wind-generated waves in coastal regions and inland waters.¹¹ It solves the evolution equation of wave action density N(x→,t;σ,θ) in space x→ and time t over frequency σ and wave direction θ,

where c→g is the group velocity, U→ is the mean current. The quantities c_σ and c_θ are the propagation velocities in spectral space (σ, θ), and S_tot is the source term. In (1), the left-hand side describes the wave kinematics, and the right-hand side accounts for various wave energy sources and sinks, including wave generation by wind, wave decay due to whitecapping, bottom friction, and depth-induced wave breaking, as well as energy redistribution through nonlinear wave-wave interactions. Recognized as a reliable coastal community wave model, SWAN has been widely applied for wave hindcasting and forecasting (Rogers et al., 2003, 2007) for the Great Lakes (Anderson et al., 2015; Mao et al., 2016; Huang et al., 2021) and coastal oceans (Dietrich et al., 2011; Huang et al., 2013; Cheng et al., 2015; Xie et al., 2016; Niroomandi et al., 2018; Xie et al., 2019).

A curvilinear grid mesh was applied in the SWAN simulation, which consists of 181 × 496 grid cells with a grid resolution of ∼1 km (Figure 1B). The spectral domain was discretized into 12 directions with 30-degree intervals and 32 frequency bands from 0.0521 to 1.0 Hz. The default values of the parameters in SWAN model were employed. The model was driven by hourly surface wind forcing from the CFSR and CFSv2.

To address the effects of ice coverage on wind-generated waves, a common technique is to mask the mesh cell with land when the local ice coverage exceeds a threshold value (e.g., Hubertz et al., 1991; Bennington et al., 2010; Tuomi et al., 2011; Anderson et al., 2015). This was realized through bathymetry modifications by changing the ice-covered areas to land points (Tuomi et al., 2019). The temporal variations of daily ice cover were from the NOAA Great Lakes Ice Cover Database described in section “Data”. Following the study in Great Lakes by Anderson et al. (2015), the value of 30% was selected as the threshold value used in this study, i.e., the wave energy was assumed to be totally dissipated where the ice coverage exceeds 30% at a given model grid cell.

A 42-year long hindcasting was performed from 1979 to 2020. The model results were compared with the measured significant wave heights (H_s) and peak wave period (T_p) at 18 buoy stations in Lake Michigan (Figure 1B). Among these buoys, S45002 (45.344N, 86.411W) was installed in the northern basin in 1979, while S45007 (42.674N, 87.026W) in the southern basin in 1981. The other 16 buoys were deployed in shallower water along the coast within less than 7 km to the shoreline in the last 10 years (Figure 1B and Table 1). Since the buoys were removed during cold seasons (usually from December to April), wave measurements were not available for this period. The buoy data availability is summarized in Table 1 with more details for the two long-term buoy records at S45002 and S45007 shown in Figure 2.

The performance of the model was assessed using the following four statistical error measures: Bias, Root Mean Square Error (RMSE), Scatter Index (SI), and R². The bias is the difference between the calculated and measured mean values and given by

where Y¯=∑i=1NYi/N is the mean of the calculated values Y_i (i = 1, 2..N) from the model and X¯=∑i=1NXi/N is the mean of the measured values X_i (i = 1, 2..N) from observation with sample size of N. The RMSE is calculated using the following formula,

The SI is defined as the normalized RMSE by the measured mean value and given by Lin et al. (2002)

The R² is the coefficient of determination defined as

The comparisons between the calculated and measured H_s and T_p were conducted for all the available buoy data from 1979 to 2020. For the purpose of demonstration, direct model-data comparisons of H_s and T_p based on hourly results for 2020 at the two offshore buoys and 12 coastal buoys (no data available from the other four coastal buoys) are shown in Figures 3, 4, respectively. Such comparisons for other years are presented in the (Supplementary Figures 1–8). The model showed consistent performance over the entire simulation period. As shown in Figure 3, the hourly wave heights oscillated in a range of 0–3 m during the year, with stronger waves generally after mid-September. Substantial spatial variability was shown across the lake. The largest wave heights were observed offshore and exceeded 2.8 m at S45002 and S45007. Coastal waves were generally weaker, yet strong waves with H_s greater than 2 m occurred on the east coast at S45024. In other coastal regions, H_s was generally below 2 m, and the H_s along the southwest coasts at S45186 was particularly small with H_s less than 1m in most of the year.

The model captured the spatiotemporal pattern of the wave characteristics (Figure 3) very well. Particularly in offshore regions, the modeled H_s agreed well with the measurements at S45002 and S45007 with R² = 0.91 and R² = 0.83, respectively. For coastal waves, the model also showed a close agreement with observation with R² values between 0.73 and 0.86 for all but two buoys (S45175 and S45186). The results showed the model explained at least 70% of the variability of the observed data around its mean. In addition, the model showed small RMSE (0.15–0.29 m) and bias (± 0.09 m) except for S45174 and S45170 with bias of –0.15 m and –0.18 m, respectively, which is likely to be caused by the underestimation for large waves. Comparisons of T_p were generally less desirable but still showed a reasonably good agreement (R² = 0.76 and R² = 0.66) at the two long-term record offshore buoys S45002 and S45007 (Figures 4A,B). At coastal buoys, the model-data comparisons of T_p were reasonably well, except for the comparison at S45175 (Figure 4I). The RMSE for T_p comparison was less than 1.08 s except for S45175 (RMSE = 1.23 s), where available observation data were very limited. Additionally, the results also showed small biases of less than 0.30 s except for S45161 and S45174, which had larger biases of 0.78 s and 0.46 s, respectively.

Long-term wave hindcasting for annual mean H_s, 90th percentile H_s, and 99th percentile H_s from 1979 to 2020 was evaluated (Figure 5). For appropriate comparisons, the annual statistics for buoy measurements and model results were calculated for the same time period based on the buoy data availability (Figure 2). This is referred to as Type M statistics in Tuomi et al. (2011), i.e., statistics calculated only from the time coinciding with the measurements. Results have shown that the model accurately captured the interannual variability of significant wave height for its mean and 90th percentile values. For example, in the northern basin (S45002), both mean and 90th percentile wave heights showed a continuous increase during 1980–1985, followed by a relatively large decrease in wave heights in 1986 before the wave heights increased again until 1990 and decreased throughout 1995. In addition, the 90th percentile H_s in the northern basin in 2007 was noticeably higher than previous and following years, which has also been accurately captured by the model. The model generally underestimated the extremely high waves (i.e., 99th percentile H_s) before 2010, particularly in the years of 1995, 2001, and 2007. Nonetheless, the model performed well in capturing the 99th percentile wave heights after 2010. This is likely due to the improved wind forcing data from the CFSv2 available from 2011.

In addition to the statistics of annually averaged wave characteristics, the scatter plots of the hourly H_s were also assessed across the lake (Figure 6, Type M statistics). Six representative sampling locations were selected, including two offshore buoys: northern basin (S45002) and southern basin (S45007), four coastal buoys that represent northeast coast (S45022), east coast (S45024), southeast coast (S45026), and southwest coast (S45186). The model explained the variability of data well, indicated by a high R² value (0.66–0.82). The SI ranged between 0.3 and 0.5, which is also similar to the SI ranges reported in other wave modeling evaluations (Lin et al., 2002; Anderson et al., 2015; Vieira et al., 2020). More importantly, the Quantile-Quantile (Q-Q) plots (magenta crosses) show that wave height distributions from model and observations were highly consistent in the offshore basins, northeast, and east coasts as Q–Q plot approximately lied close to the line y = x. In the southeast and southwest coasts, the model was able to accurately reproduce weak wave conditions (i.e., H_s < 1 m) but underestimated large waves (e.g., 99th percentile H_s), particularly at S45026.

There are several ways to calculate wave statistics in the presence of a seasonal ice cover as laid out in detail by Tuomi et al. (2011). For example, we have used Type M statistics in section “Model Configuration” for model validation. Other main types of statistics described in Tuomi et al. (2011, 2019) include Type I (statistics with ice time included when H_s = 0 m), Type F (ice time excluded completely from the calculations), Type N (hypothetical no-ice conditions). In the following sections, we use the type I analysis for H_s unless noted otherwise. Namely, wave analysis performed at any given location utilizes the entire 42-year simulation results including the ice time with H_s = 0 when the local ice cover exceeds the threshold value of 30%. In such a way, it allows us to directly assess the impact of ice cover change on the wave climate.

To understand the seasonality and spatial variability of wave height, the modeled 42-year (1979–2020) mean lake-wide average significant wave height H_s, 90th percentile, and 99th percentile for each month, respectively, are provided in Figure 7. It is clear that the wave height variability closely aligned with the wind speed at 10 m height (U₁₀) except in the mid ice season. Winds were the strongest in January but the presence of ice cover reduced the wave height compared to December. In February, decreasing U₁₀ and increasing ice cover both helped reduce waves. However, in March and April, decreasing U₁₀ and decreasing ice cover offset their negative and positive impacts on wave growth, leading to an insignificant change in H_s. With the disappearance of ice after April, the wave height dropped rapidly from May to August in response to the decreasing wind during summertime. The increasing wind in the fall enhanced wave development, and the wave height reached the annual peak in December when the winds became strong in wintertime and ice had yet to develop.

Since the buoy data were usually not available from December to April, the observation data provided an incomplete picture of the wave characteristics and would underestimate the annual wave conditions. The model results are, therefore, essential to provide a comprehensive understanding of the wave climate in the entire lake and throughout all seasons. This was depicted by comparing the wave patterns averaged over the ice-free season, the ice season, and the entire year (Figure 8). We note that the ice season in this work is defined as December to April, which is based on the monthly climatology of ice cover in Figure 7. It is possible that ice may not be present at some location during ice season or may be present during the ice-free season over the 42-year simulation time. This should not be confused with the time of ice cover at a given specific location. All results in Figure 8 are based on the Type I analysis. Waves were stronger offshore than along coastal regions. On a close look, meridional differences in mean H_s between the northern and southern basins were noticeable during the ice season, with stronger waves in the southern basin resulting from the dominant northwesterly wind during the winter (Figure 8A2). On the other hand, the zonal difference in large waves (i.e., 99th percentile H_s) between west and east became more evident, showing stronger waves in the eastern part of the lake (Figures 8C1–C3).

To address whether linear trends of wave climate existed across the lake that may facilitate the wave climate projection, we investigated the spatial pattern of the trends of H_s for the entirety of Lake Michigan. Linear regression analyses of H_s were conducted regarding the mean, 90th percentile, and 99th percentile of annual averaged H_s for the past 42 years (1979–2020). The MATLAB function Regstats was used to obtain all the regression coefficients based on the least square estimator and their 95% confidence intervals. Regions with a statistically significant trend at 95% confidence level are marked in Figure 9. The mean H_s and 90th percentile H_s showed a statistically significant increase in some small areas in the northeast, east, and southeast of the lake. For the 99th percentile H_s, a statistically significant increase existed only in smaller regions in the northeastern and eastern coasts. In contrast, a downward trend was observed within a band from the west coast to the south basin. The results have shown that there was no simple linear trend of long-term H_s change for the majority (> 90%) of the lake. Similarly, we applied the same linear trend analysis to H_s of successive January(s), February(s), …, December(s) of the 42 years, respectively. The results were consistent with Figure 9 that no linear trend could be established for most areas of the lake.

The above analysis suggests the inadequacy of linear trend analysis used in previous studies (Olsen, 2019; Jabbari et al., 2021). It also reinforces the complexity of wave climate in Lake Michigan. A 10-year trailing moving mean was therefore applied to H_s to remove seasonal and annual variability, focusing on long-term wave patterns (Figure 10). Simulated waves at all buoy locations suggest there was no clear trend and instead, regime shifts were identified. Specifically, there were trends of increasing H_s around 2010–2020 along with rising water levels at 14 out of the 16 representative buoy stations, except S45013 and S45186 at the southwest coast, where the increasing trends after 2010 were less certain. However, it should be noted that S45187 and S45174 were also located on the southwest coast and showed clear trends of increasing wave height. Furthermore, all stations showed low H_s around the period of 2000–2010 (coincident with low lake levels). In addition, 14 stations (except S45020, S45022) illustrated a downward trend around 1990–2000. These patterns resembled the annual water level change in the region with the correlation coefficients greater than 0.8 at five stations, between 0.7 and 0.8 at two stations, and between 0.6 and 0.7 at four stations (Figure 10). Increasing H_s possibly reflected increases in storm intensity and frequency, preceding the rising lake-wide water level.

Water level changes do not directly influence wave height significantly, except very shallow nearshore regions where the water level change would impact the wave shoaling and breaking. However, results in Figure 10 have revealed the correlation between changes in water level and wave height across the lakes in both offshore and coastal regions that are further away from nearshore zones. The correlation between water level and wave height does not suggest any causal relationship between the two variables. Instead, such correlation suggests they were both influenced by atmospheric and regional climate patterns. On one hand, the wave height is known to be primarily controlled by local wind speed, duration, and fetch in the Great Lakes (Pore, 1979; Donelan et al., 1985), while in coastal oceans swells also play an important role (Aboobacker et al., 2011; Sabique et al., 2012; Wandres et al., 2020). The H_s and wind speed at a 10 m height showed a strong correlation with a correlation coefficient of r≥0.90 at 12 out of the 16 stations (not shown). On the other hand, the primary drivers of water levels in the Great Lakes are runoff, over-lake precipitation, connecting channel flow, and evaporation (Deacu et al., 2012). Regional climate assessment showed that more intense storms due to hydrologic intensification have been observed in the Great Lakes region in the past decades. The amount of precipitation falling in the heaviest 1% of storms increased by 35% in the Great Lakes region from 1951 through 2017 (USGCRP, 2017). Therefore, the changing pattern in storms, which influences waves (instantly through increased wind stress on the surface) and water levels (following, with a lag through precipitation and runoff), is likely to be a driving factor responsible for the observed correlation between the wave height and water level.

With a changing regional climate, the Lake Michigan water level has shown three regimes in the past four decades with a high water level between 1980 and 1998, remained an extremely low water level for more than a decade during 2000–2013, and followed by a drastic and rapid increase (water level increased by 2 m within 6 years) during 2014–2020 (Figure 11). The shifts among these three regimes occurred around 1998 and 2014, corresponding to two strong El Niño events, i.e., the 1997–98 and 2014–16 events. To further validate the hypothesis on coherent patterns of the waves, water level, and their possible connection to the regional climate change including storminess and ice reduction; we analyzed the differences in large waves (99th percentile H_s) during these three regimes defined by water level change, as they are a more suitable indicator of storm impacts than mean waves height. The changes in 99th percentile H_s from the first regime (1979–1998) to the second regime (1999–2013) and from the second regime (1999–2013) to the third regime (2014–2020) are shown in Figure 12.

From the period of 1979–1998 to the period of 1999–2013, the wave height (99th percentile H_s) averaged over each period, decreased by 0.1–0.3 m in the majority of the lake (Figure 12A1). This is consistent with the shift from high water to low water levels (Figure 11). The exception is in the northernmost region where the wave height increased by < 0.1 m due to decreasing ice cover, which is discussed in the following section. During the ice season, a similar pattern of wave change was observed (Figure 12A2), with changes in the magnitude of H_s becoming more drastic. Unlike the pattern during ice season, the changing pattern during the ice-free season showed a somewhat opposite pattern with an increase in wave height in the eastern lake and southern lake while a decrease was observed in the western and northern lake (Figure 12A3). These results suggest the complexity of spatiotemporal patterns of storm changes, indicating the possibility of declined winter storms and increased summer storms during this period. However, the annual wave pattern is significantly influenced by the wave pattern during the ice season, highlighting the essential role of winter wave conditions in wave climate.

In contrast, from 1999–2013 to 2014–2020, wave heights increased in most regions of the lake by 0.1–0.3 m, including the entire northern lake and the eastern part of the southern lake (Figure 12B1). In addition, the increases in wave heights during the period of 2014–2020 were consistent across the ice season, ice-free season, and annual average, showing an apparent increase in wave heights that were well aligned with rapid water level rise. This reinforces the hypothesis of coherent patterns of the waves, water level, and their connection to more frequent and intensified storms observed in the Great Lakes in the past.

In this section, we made an attempt to examine if the changes in wave climate can be detected in wave height extreme value and its impact on design wave heights corresponding to different return periods. Independent events were selected from the continuous time series of H_s using the peaks-over-threshold (POT) method. To ensure the independence of events, the selected events were recommended at least 48 h apart (Caires and Sterl, 2005; Aarnes et al., 2012; Anderson et al., 2015; Velioglu Sogut et al., 2018). Following Méndez et al. (2006) and Niroomandi et al. (2018), a minimal time interval of 72 h between events was adopted in this study. The threshold value for individual wave events was set as 3 m, which is comparable to the 99th percentile for H_s in the offshore station 45002 during the 42-year simulation. Sensitivity analysis was also conducted by selecting a threshold of 3.3 m (99.5th percentile), results remained nearly the same. The return period analysis was conducted with the selected wave events by fitting the Generalized Pareto distribution (GPD). The cumulative distribution function for GPD was given by Coles (2001),

where y values are threshold excesses, σ is the scale parameter, and ξ is the shape parameter. The parameters σ and ξ were obtained by maximum likelihood estimates. Then, the return level was estimated by Coles (2001),

where μ is the threshold, n_y is the average number of the POT selected events per year, and N is the return period in years.

Design wave heights at the offshore buoy location S45002 corresponding to different return periods were calculated repeatedly using simulation data from three periods (Figure 13), namely, Period I: 1979–1998 (regime 1); period II: 1979–1998 plus 1999–2013 (regime 1 plus regime 2), and period III: 1979–1998 plus 2014–2020 (regime 1 plus regime 3). The 100-year return level of H_s during period I was estimated to be 6.05 m. Using the data in period II, the 100-year return level of H_s reduced slightly to 5.87 m, which reflects the impact of the reduced H_s during regime 2 as discussed before. Similarly, with the data in period III, the 100-year return level of H_s was elevated to 5.98 m. This again reflects the increased H_s during regime 3. Consistently, the events exceeding the threshold decreased during period II and increased during period III (Figure 13). However, these estimated 100-year return levels only vary moderately (<2% variation relative to the mean of these values = 5.96 m) using data from the three periods, suggesting the robustness of the estimates. It also suggests that the changes in wave climate during the past four decades have had no significant impact on design wave weight and extreme value analysis.

Ice cover can reduce wind stress and resulting wave energy. In addition, waves propagating in an ice field are reduced by a combination of scattering and dissipation. Due to the orientation of Lake Michigan, the northernmost part of the lake is characterized by a much higher ice cover percentage and longer ice duration, while the southern region has much less ice cover during the winter. Therefore, the trend and variability of ice cover and its impact on waves are much more significant in the northern lake. Recall that the northernmost part of the lake showed an increase in large (e.g., 99th percentile) wave heights in ice season (Figure 12A2) from the period of 1979–1998 (high water level) to the period of 1999–2013 (low water level), which is opposite to the downward trend observed in the most areas of the lake which had no or low ice cover. Such a difference was caused by the significant reduction in ice cover in favor of wave development that occurred in the northern lake during the period of 1999–2013 in comparison to the period of 1979–1998 (Figures 14A,B). It also explains why the northernmost part of the lake could show a statistically significant increase in wave height over the past 42 years, even with linear regress analysis. Lastly, results (Figures 12B2, 14C) suggest that the increased winter storms accompanied by the rising water level during the period of 2014–2020 were the dominant forcing for the increase in large waves, which overshadowed the negative influences of the moderate increase in ice cover in the northern lake on wave development.

Under projected regional warming, the observed reduction in ice coverage and ice duration in the Great Lakes are expected to continue into the future (Wang et al., 2012; Mason et al., 2016; USGCRP, 2017). To demonstrate the influence of changing ice conditions on the wave climate in Lake Michigan, we selected two representative historical years, which have low and high ice covers (2013, 2014) for case studies with both Type I and Type N analysis (Figure 15). The wind speeds averaged over the lake during the ice seasons were similar in the two years, with U₁₀ = 7.78 m/s in 2013 and U₁₀ 7.50 m/s in 2014. However, the averaged ice coverage in 2013 was only 11% compared to a much higher percentage of ice cover of 51% in the cold year of 2014. While the wind speed in 2014 ice season was only slightly lower than that in 2013 during the ice season (Figures 15A2,B2), the mean H_s in 2014 ice season was 0.2–0.8 m, which was significantly lower than that in 2013 during which the majority of the lake showed mean H_s between 0.8 and 1.2 m (Figures 15A3,B3, Type I statistics). In a warming scenario of no ice (the wave simulations in these 2 years were re-run under the condition without ice cover in the lake to demonstrate the impact of disappearing ice on the wave climate, Type N statistics), mean H_s would increase noticeably in the northernmost part of the lake, a region influenced the most by the ice cover, in the simulation of 2013. In the simulation of the cold year of 2014, a substantial increase by 50% in mean H_s would occur in the majority of the lake, with the most significant increase in the northern lake, where ice cover was the largest (Figures 15A4,B4). This again illustrated the influences of decreasing ice on wave development.