2. Methods
To characterize changes in lake level variability, we analyzed lake level data compiled by the Great Lakes Coordinating Committee (Coordinating Committee on Great Lakes Basic Hydraulic Hydrologic Data, n.d.). Monthly mean lake level data are calculated by combining observations from a suite of gauges around the Great Lakes basin maintained by the U.S. National Oceanic and Atmospheric Administration (NOAA) and Canadian Hydrographic Service (CHS) (Fry et al., 2022). These data span from 1918 to the present. In addition to lake level data, our analysis uses Great Lakes monthly hydrologic data curated by the NOAA Great Lakes Environmental Research Laboratory (GLERL) for net basin supply variables: overlake precipitation, overlake evaporation, and runoff.
The NOAA-GLERL monthly hydrometeorological database (GLM-HMD) is the first dataset to assimilate hydrometeorological measurements across both the U.S. and Canadian portions of the Great Lakes basin dating back to the early and mid-1900s (Hunter et al., 2015). The GLM-HMD uses a Thiessen weighting approach to compute overlake precipitation from daily monitoring stations (Hunter et al., 2015), including stations up to 50 km from the lake basin depending on gauge density, though stations over the lake surface or on islands are sparse. It should also be noted that station availability expanded rapidly from the early 1900s before peaking from the late 1940s to the present (see fig. 2, Hunter et al., 2015) and GLM-HMD overlake precipitation calculations are thus available beginning in 1940. Overlake evaporation is computed by the Large Lake Thermodynamics Model (LLTM) as described in Croley (1989), and is available from 1950 to present. The LLTM is a 1-dimensional thermodynamics model that computes overlake evaporation by simulating the energy balance above the lake surface, vertically-summed heat storage through the lake columns, and aerodynamic evaporation (Fry et al., 2022). GLM-HMD estimates of runoff are calculated from streamflow measurements from the United States Geological Survey and Water Survey of Canada gauges based on subbasin area with extrapolations for unrepresented subbasin area and ungauged subbasins (Hunter et al., 2015). Fry et al. (2013) determine that this approach provides reliable estimates of total discharge to the lakes. In preparation to examine trends in interannual variability, we transformed monthly data to annual data by taking an annual mean (lake levels) or annual sum (precipitation, evaporation, runoff) on the component time series as appropriate. Finally, we note that we considered alternative data sources, notably output from the Large Lake Statistical Water Balance Model (L2SWBM) (Gronewold et al., 2020), and a comparison of annual values from GLM-HMD and L2SWBM produced correlations of at least 0.97 for all lake-component combinations.
Here, we define interannual variability as the standard deviation in a time series over a 13-year moving window. We select a 13-year window width for two reasons. First, this width is longer than the cycles of known 8- and 12-year modes of lake level variability (Hanrahan et al., 2014; Cheng et al., 2021). Interannual lake level fluctuations are driven both by modes of climate variability at various periodicities, with notable periodicities of 8, 12, and 30 years (Hanrahan et al., 2014; Cheng et al., 2021), as well as by long-term trends in net basin supply components. Relevant modes of climate variability include the Pacific Decadal Oscillation (PDO), Atlantic Multidecadal Oscillation (AMO), North Atlantic Oscillation (NAO), El Nino-Southern Oscillation (ENSO), and the Pacific-North American Oscillation (PNA) (Wang et al., 2018; Ghanbari and Bravo, 2008). Second, we examined the dependence of trends through a sensitivity analysis using window widths ranging from 5 to 25 years. This revealed a consistency in trends in interannual variability for window widths greater than 12 years despite inconsistencies in trends across shorter window widths. We focus our analysis on the period from 1970 to the present. While this truncates the length of data available for time series analysis, this is a period of increasing influence of anthropogenic climate change. Notably, there is evidence of an emerging intensification of the hydrologic cycle early in this five-decade span or just prior to it (IPCC, 2021), an intensification which is expected to continue (Seager et al., 2014). However, we include the full extent of available data in Figure 1 to provide historical context for recent trends.
To calculate the time dependence of lake levels and NBS components, we use two non-parametric measures: the Kendall rank correlation coefficient and the Trend-Free Pre-Whitened Mann-Kendall trend test. We selected these tests due to the appropriateness of their application on time series with high levels of autocorrelation; this restriction ruled out more traditional linear regression techniques (e.g., the Ordinary Least Squares and conventional Mann-Kendall trend tests) whose statistical assumptions are not satisfied by the time series under analysis. The Kendall rank correlation coefficient is a rank-based test that determines if two variables are related regardless of distributions. The Trend-Free Pre-Whitened Mann-Kendall trend test is a modified version of the conventional Mann-Kendall trend test designed to overcome inaccuracies in analysis due to autocorrelation (Yue et al., 2002); both the modified and conventional Mann-Kendall trend tests are valid independent of the distribution of the underlying data. For ease of reading, we will refer to the Trend-Free Pre-Whitened Mann-Kendall trend test simply as the Mann-Kendall trend test throughout this study. In summary, to assess the robustness in trends of interannual variability, we apply a Mann-Kendall trend test to the 13-year moving window of standard deviation across lake levels and NBS components over the past fifty years.