Wave Climate Associated With Changing Water Level and Ice Cover in Lake Michigan

Detailed knowledge of wave climate change is essential for understanding coastal geomorphological processes, ecosystem resilience, the design of offshore and coastal engineering structures and aquaculture systems. In Lake Michigan, the in-situ wave observations suitable for long-term analysis are limited to two offshore MetOcean buoys. Since this distribution is inadequate to fully represent spatial patterns of wave climate across the lake, a series of high-resolution SWAN model simulations were performed for the analysis of long-term wave climate change for the entirety of Lake Michigan from 1979 to 2020. Model results were validated against observations from two offshore buoys and 16 coastal buoys. Linear regression analysis of significant wave height (Hs) (mean, 90th percentile, and 99th percentile) across the entire lake using this 42-year simulation suggests that there is no simple linear trend of long-term changes of Hs for the majority (>90%) of the lake. To address the inadequacy of linear trend analysis used in previous studies, a 10-year trailing moving mean was applied to the Hs statistics to remove seasonal and annual variability, focusing on identifying long-term wave climate change. Model results reveal the regime shifts of Hs that correspond to long-term lake water level changes. Specifically, downward trends of Hs were found in the decade of 1990–2000; low Hs during 2000–2010 coincident with low lake levels; and upward trends of Hs were found during 2010–2020 along with rising water levels. The coherent pattern between the wave climate and the water level was hypothesized to result from changing storm frequency and intensity crossing the lake basin, which influences both waves (instantly through increased wind stress on the surface) and water levels (following, with a lag through precipitation and runoff). Hence, recent water level increases and wave growth were likely associated with increased storminess observed in the Great Lakes. With regional warming, the decrease in ice cover in Lake Michigan (particularly in the northernmost region of the lake) favored the wave growth in the winter due to increased surface wind stress, wind fetch, and wave transmission. Model simulations suggest that the basin-wide Hs can increase significantly during the winter season with projected regional warming and associated decreases in winter ice cover. The recent increases in wave height and water level, along with warming climate and ice reduction, may yield increasing coastal damages such as accelerating coastal erosion.

Detailed knowledge of wave climate change is essential for understanding coastal geomorphological processes, ecosystem resilience, the design of offshore and coastal engineering structures and aquaculture systems. In Lake Michigan, the in-situ wave observations suitable for long-term analysis are limited to two offshore MetOcean buoys. Since this distribution is inadequate to fully represent spatial patterns of wave climate across the lake, a series of high-resolution SWAN model simulations were performed for the analysis of long-term wave climate change for the entirety of Lake Michigan from 1979 to 2020. Model results were validated against observations from two offshore buoys and 16 coastal buoys. Linear regression analysis of significant wave height (H s ) (mean, 90th percentile, and 99th percentile) across the entire lake using this 42year simulation suggests that there is no simple linear trend of long-term changes of H s for the majority (>90%) of the lake. To address the inadequacy of linear trend analysis used in previous studies, a 10-year trailing moving mean was applied to the H s statistics to remove seasonal and annual variability, focusing on identifying long-term wave climate change. Model results reveal the regime shifts of H s that correspond to long-term lake water level changes. Specifically, downward trends of H s were found in the decade of 1990-2000; low H s during 2000-2010 coincident with low lake levels; and upward trends of H s were found during 2010-2020 along with rising water levels. The coherent pattern between the wave climate and the water level was hypothesized to result from changing storm frequency and intensity crossing the lake basin, which influences both waves (instantly through increased wind stress on the surface) and water levels (following, with a lag through precipitation and runoff). Hence, recent water level increases and wave growth were likely associated with increased storminess observed in the Great Lakes. With regional warming, the decrease in ice cover in Lake Michigan (particularly in the northernmost region of the lake) favored the wave growth in the winter due to increased surface wind stress, wind fetch, and wave transmission. Model INTRODUCTION Lake Michigan is the second largest of the Laurentian Great Lakes by water volume (about 4,900 km 3 ) and the third largest by surface area (about 58,000 km 2 ) with a maximum length (north-east) of about 494 km and a maximum width (east-west) of approximately 190 km ( Figure 1A). Coastal communities of Lake Michigan are confronted with challenges to manage the 1,640 miles shoreline under changing coastal environments. Wind waves have a significant influence on Lake Michigan's geomorphological processes by inducing lakebed sediment resuspension (Schwab et al., 2006), modifying outer and inner nearshore morphologies (Mao et al., 2016), and impacting bluff recession and shoreline erosion (Meadows et al., 1997;Swenson et al., 2006;Lin and Wu, 2014;Anderson et al., 2015). Large waves have also caused coastal wetland losses, yielding the reduction of habitats and species diversity (Thomasen et al., 2013). Understanding the wave climate changes in Lake Michigan is essential for coastal protection and restoration to address the challenges of excessive shore and beach erosion, severe bluff and dune recession, and the damages to coastal protection structures by climate or weather hazards (Panchang et al., 2013). Within the past decade, Lake Michigan's water level has rapidly transitioned from decadal-long record low levels to record highs (Gronewold and Rood, 2019) due to the hydrologic intensification. This may have potentially exacerbated wave conditions that may cause more severe damage to beaches, shoreline, coastal infrastructure, communities, and ecosystems. Therefore, there is an urgent need for assessing the present change in wave climate and its future trend.
While numerous studies have been performed to investigate the long-term wave statistics in oceanic and coastal regions (e.g., Tuomi et al., 2011Tuomi et al., , 2019Camus et al., 2014;Bromirski and Cayan, 2015;Cheng et al., 2015;Wang et al., 2015;Niroomandi et al., 2018), there are only a few systematic lake-wide wave climate studies for the Great Lakes (e.g., Hubertz et al., 1991;Jensen et al., 2012), partly due to the limited long-term wave records. In Lake Michigan, all coastal buoys were deployed during 2010-2018 (Table 1). Long-term wave data with more than 30-year records, suitable for wave climate analysis (Trewin, 2007), are only available at two offshore National Data Buoy Center (NDBC) buoys [Station 45002 (hereafter S45002) and S45007, data available since 1979 and 1981, respectively], in the northern and southern basins ( Figure 1B). In addition, the buoy data are not available during cold seasons (usually from December to April, Figure 2). Based on the data at the two offshore buoy stations, Olsen (2019) found that the measured yearly mean significant wave heights (H s > 0.25 m from May to October) from 1981 to 2017 were decreasing at the rates of 1.5 mm/yr and 3.8 mm/yr at S45002 and S45007, respectively, with 95% confidence intervals using Theil-Sen regressions. Recently, Jabbari et al. (2021) investigated the trend of H s in the months of August  at S45007 using linear regression and found an increasing trend of 2.4 mm/yr. In their data analysis, the missing wave data were filled using the empirical relation between H s and the measured wind speed (Carter, 1982). The contradiction between Olsen (2019) and Jabbari et al. (2021) indicates that the trend of H s in Lake Michigan may be nonlinear or even more complex. Additionally, the limited number of buoys (two for long-term analysis) is insufficient for accurate analysis of wave climate over the entire lake, especially under severe storm conditions (Mao et al., 2016).
To fill the gaps of buoy data and develop the spatial distribution of wave climate, numerical models were often used (Caires and Sterl, 2005;Jensen et al., 2012;Bromirski and Cayan, 2015). Using the spectral wave model WAM Cycle 4.5.1C, the U.S. Army Corps of Engineers carried out long-term wave hindcasts for all the five Great Lakes under the Wave Information Studies (WIS) project (Hubertz et al., 1991;Jensen et al., 2012). For Lake Michigan, the WIS provided wave hindcasts from 1979 to 2017 at 490 nearshore stations around the lake. With the WIS data at nearshore stations, Olsen (2019) found there were only a few stations in the southwest and two stations in the southeast of the lake that had statistically significant decreasing trends for the yearly mean of H s > 0.25 m while there were no significant trends for the majority of WIS stations. On the other hand, using the 32-year DWAVE model (Resio and Perrie, 1989), an earlier model version used in the WIS project from 1956 to 1987, Meadows et al. (1997) established a correlation between the Great Lakes water levels and incident wave energies, having hypothesized that they were associated with large-scale shifts in atmospheric storm tracks. This again indicates the complex changes in wave climate influenced by regional climate change.
Since the study of Meadows et al. (1997), the water level of Lake Michigan in the past four decades has shown strong interannual variability. The lake level declined from a high water level in 1987 to a level below the long-term mean in 1999 and remained very low (setting an all-time low in January 2013) for more than a decade until 2014. Since then, the lake level has increased rapidly by ∼2 meters in 6 years. In 2020, the water level continually broke the monthly record high through January to August. Meanwhile, Lake Michigan has experienced some of the biggest waves with wave height reaching 6.5 m during fall storm events in the past decade (e.g., Sep 30, 2011, 1 October 31, 2014, 2 FIGURE 1 | (A) Bathymetry and location of Lake Michigan in the United States. (B) Computational domain and curvilinear grid used for the wave hindcasting. The triangles indicate the locations of MetoOcean buoys S45002 and S45007 with long-term data from 1979 to 2020 and 1981 to 2020, respectively, while the circles indicate the locations of coastal buoys with short-term (later than 2010) data availability. Figures 1, 8 , 9, 12, 14, 15 are produced following the guidance for colormap selection by Thyng et al. (2016).
March 6, 2020 3 ). Additionally, the ice coverage during winter has significantly reduced in the Great Lakes (Assel et al., 2003;Wang et al., 2012) in response to regional warming trends. As ice coverage affects the wind waves by reducing wind stress and dissipating wave energy (Bennetts et al., 2010;Doble et al., 2015;Sutherland and Rabault, 2016;Kohout et al., 2017;Voermans et al., 2019;Yiew et al., 2019;Bai et al., 2020), reduced ice coverage is expected to favor the wave growth. The combined impacts of these factors on the wave climate in Lake Michigan have not been fully understood.
The objective of this study was to investigate the response of long-term wave climate for the entirety of Lake Michigan to the changes in regional climate factors using a numerical wave model. The model was validated against the measured wave data from offshore and coastal MetOcean buoys. The validated model was then applied to investigate the spatial and seasonal characteristics 3 https://www.weather.gov/lot/2020March6_lakeshoreflooding of the wave climate in Lake Michigan with a focus on wave climate change. To investigate the long-term changes in H s , the validity of conventional linear regression used in previous studies was examined. Furthermore, after removing the seasonal and annual variability with the 10-year moving average, we identified the regime shifts of H s and its correlations with wind speed, water level, and ice cover. Finally, case studies were conducted to explore the impact of decreasing ice cover on the wave climate.

Data
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) 4 dataset with a spatial grid resolution of ∼0.3 degree for 1979-2010 and the upgraded Climate Forecast System Version 2 (CFSv2) 5 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(Saha et al., , 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(Saha et al., , 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, 6 which includes the Great Lakes Ice Atlas 7 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. 8 These lakewide 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. 9 These buoys consist of a 1.1 m-diameter conical hull with a six-degree-offreedom 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, 10 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.

Model Configuration
Prevalent numerical models for wave hindcast are the thirdgeneration 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. 11 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(Rogers et al., , 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.

Model Validation
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 2 . The bias is the difference between the calculated and measured mean values and given by where Y = N i = 1 Y i /N is the mean of the calculated values Y i (i = 1, 2..N) from the model and X = N i = 1 X i /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 2 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 2 = 0.91 and R 2 = 0.83, respectively. For coastal waves, the model also showed a close agreement with observation with R 2 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 2 = 0.76 and R 2 = 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 2 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.

Wave Climate in Lake Michigan
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. (2011Tuomi et al. ( , 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 42year 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 10 ) 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 10 and increasing ice cover both helped reduce waves. However, in March and April, decreasing U 10 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).

Wave Regime Shift
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 . 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.
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.

Impact on Wave Extreme Value
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: 1979II: -1998II: plus 1999II: -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.

Effects of Ice Reduction on Wave Climate
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 10 = 7.78 m/s in 2013 and U 10 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.

CONCLUSION
In this study, we investigated the wave climate associated with changing water level and ice cover in Lake Michigan. The main findings include: (1) long-term wave changes were poorly explained by a simple linear trend. Rather, we identified the several regime shifts of H s that correspond to long-term lake water level changes. The coherent patterns of the waves and water levels were hypothesized to be due to the changes in the storms crossing the entire lake basin, which influenced the changes in waves (instantly through increased wind stress on the surface) and water levels (following, with a lag through precipitation and runoff). (2) The seasonal ice cover has a decisive impact on the wave climate. In our case studies, simulations with hypothetical no-ice conditions reveal that loss of ice cover can increase basin-wide mean wave height by 50% during the winter season. Under the projected climate, the regional air temperature is projected to increase by 3.3 • C-6.0 • C by the end of the century, depending on warming scenarios (Zhang et al., 2020). Correspondingly, the ice cover in the Great Lakes is expected to continue to decrease. Lake Michigan's recently increasing wave heights and water levels, along with a warming climate and ice reduction, may yield increasing coastal damages such as accelerating coastal erosion and enhanced sediment transport.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s. a grant (NA18NOS4190162) from the National Oceanic and Atmospheric Administration, U.S. Department of Commerce.

ACKNOWLEDGMENTS
This is the contribution 83 of the Great Lakes Research Center at Michigan Technological University. The Michigan Tech high performance computing cluster, Superior, was used in obtaining the modeling results presented in this publication.