Changes in Climatological Variables at Stations around Lake Erie and Lake Michigan

Abstract

Climatological variables undergo changes over time, and it is important to understand such dynamic changes at global, regional, and local levels. While global and regional studies are common in the study of climate, such studies at a local level are not as common. The aim of this article is to study temporal changes in precipitation, snowfall, and temperature variables at specific stations located on the rims of Lake Erie and Lake Michigan. The identification of changes is carried out by applying change-point analysis to precipitation, snowfall, and temperature data from Buffalo, Erie, and Cleveland stations located on the rim of Lake Erie and at Chicago, Milwaukee, and Green Bay stations located on the rim of Lake Michigan. We adopt mainly the Bayesian information criterion (BIC) method to identify the number and locations of change points, and then we apply the generalized likelihood ratio statistic to test for the statistical significance of the identified change points. We follow this up by finding 95% confidence intervals for those change points that were found to be statistically significant. The results from the analysis show that there are significant changes in precipitation, snowfall, and temperature variables at all six rim stations. Changes in precipitation show consistently significant increases, whereas there is no similar consistency in snowfall increases. Temperature increases are generally quite sharp, and they occur consistently around 1985. Overall, upon combining the amounts of changes from all six stations, the average amount of change in annual average temperature is found to be 0.96 °C, the average percentage of change in precipitation is 16%, and the average percentage of change in snowfall is 17%. The changing local climatic conditions identified in the study are important for local city planners, as well as residents, so that they can be well prepared for changing climatic scenarios.

1. Introduction

The Laurentian Great Lakes, consisting of five lakes—Superior, Michigan, Huron, Erie, and Ontario—are the largest group of freshwater lakes on Earth. They contain 21% of the world’s volume of fresh surface water [1]. The Great Lakes represent a major natural resource for the United States and Canada, and they have played a vital role in the development of the industrial heartland of North America. Nearly 40 million people rely on the Great Lakes for drinking water, food, work, and recreation. The Great Lakes have a drainage area of 770,000 km² and a water surface area of 244,000 km² in the United States and Canada [2]. Lake Erie, the shallowest (average depth: ~19 m) of the Great Lakes, is home to one-third of the total human population of the Great Lakes basin [1]. These lakes have been affected by climate change in several ways, including increased surface water temperatures, longer summer stratification-related hypoxia (dissolved oxygen concentrations < 2 mg L⁻¹), and the increased occurrence of weather extremes and harmful algal blooms (HAB), which degrade the water quality [1].

A number of studies on the Great Lakes’ hydro-climatic variables have generally shown a decrease in water supplies and water levels [3,4]. Lofgren [4] analyzed the net basin supply (NBS) for Lake Superior using the coupled hydrosphere–atmosphere research model (CHARM), a regional climate model, and they predicted increases, as well as decreases, in the annual mean NBS. Yee et al. [5] reported a rapid drop in lake levels during 1987–1988. They concluded that the decrease in water levels was due to factors that included a decrease in precipitation and runoff, an increase in evaporation, and an increase in lake outflows. Mortsch et al. [6] investigated climate change impacts on the hydrology of the Great Lakes and predicted significant temperature increases in spring and winter. Ehsanzadeh et al. [2] reported that the hydro-climatology of the Great Lakes is characterized by nonstationary behavior. Based on change-point analysis, they found that precipitation and runoff were on a decreasing course, following an increasing trend in the early twentieth century. It is notable that this is the only article we found that implemented change-point methodology to study changes in the Great Lakes’ characteristics. McBean and Motiee [7] analyzed the precipitation, temperature, and streamflows of the Great Lakes for trends and found statistically significant increases in precipitation and streamflows over the period of 1930–2000. Recently, Briley et al. [8] reported that Lake Superior’s surface waters have warmed more rapidly than nearby air temperatures. Similar lake surface temperature warming has been observed worldwide [9]. Van Cleave et al. [10] reported that ice cover on the Great Lakes undergoes non-linear regime shifts and that a decline started in 1998, a year with a strong ENSO event.

Zhang et al. [11] reviewed literature that focused on climate projections derived from a variety of climate change models [12,13,14]. For example, Hayhoe et al. [12] summarized a number of climate change predictions over the Great Lakes region, focusing on specific issues and/or locations such as spring freezing and cold snaps, and the city of Chicago. Gula and Peltier [13] and d’Orgeville et al. [14] combined a land climate model with a simple lake model to generate future climate scenarios in the region.

It is also relevant to note here the impacts of changes in the Great Lakes region on the economic and climatic conditions of the land surrounding the lakes. In their recent study, Briley et al. [8] noted that large lakes can have significant impacts on regional and local climates, generating much different weather and climate conditions than when lakes were not present. In lake–atmosphere–land systems, local energy and hydrologic cycles get modified as conditions at the lake surface interact with the overlying atmosphere and nearby land surfaces. Temperature differences between the lake surface and overlying air drive lake effects such as lake breezes and enhanced lake-effect precipitation. Large lakes supply water to nearby cities, and hence, the budgets, commerce, and ecosystems of such cities are directly related to the behavior of the lake–atmosphere–land system.

The above review makes it clear that significant changes to the Great Lakes have been happening and this makes it relevant that the impacts of such changes be carried out on climatic variables in the cities surrounding these lakes. Often, such studies inform local residents about the nature of climatic changes that might be happening in their own neighborhoods, and these local changes may not be aligned with global changes. For example, the effect of Lake Erie on the snowfall in Buffalo is well known, and hence, any changes in Buffalo’s snowfall are likely to be influenced by changes that might be happening in Lake Erie [15]. A recent study based on measurements from 10,000 stations found high spatial variability in extreme precipitation [16]. The identification of high spatial variability suggests the need for more regional and local studies in order to understand the local factors that influence precipitation dynamics at a very local level.

The objective of this article is to carry out a study of changes in climatic variables at a local level. Specifically, the article focuses on studying changes in the climatic variables temperature, precipitation, and snowfall at six stations located on the rims of Lake Erie and Lake Michigan. The stations chosen in this study are part of the Automated Surface Observing Systems (ASOS) network, and the specific stations selected are as follows: Buffalo (BUF), Erie (ERI), and Cleveland (CLE), which are located on the southern rim of Lake Erie, and Chicago (MDW), Milwaukee (MKE), and Green Bay (GRB), which are located on the western rim of Lake Michigan. For each station, we considered nine specific climatic factors representing temperature, precipitation, and snowfall. Detailed information about the nine variables and the data collected at each of the stations is described in Section 2. Currently, we wish to state that this study is the first of other studies we plan to pursue on climatic changes that occur at stations on the rims of the five Great Lakes. Our final goal is to fully understand the effect of changes to the Great Lakes on climatic variables at stations on the lakes’ rims. Towards this end in this first study, we only analyzed the nature of changes at just the six stations mentioned above. Here, we do not perform a similar study on any of the lake factors, and hence, we do not attribute any direct lake effect on changes identified in this study. A larger study is being pursued as a follow-up to this initial study. In the larger study, we also plan to include changes at stations that are away from the lake rims so that we can delineate the differences in changes occurring at stations on the lake rims from those that are farther away. Of course, stations located on the rims of the other three lakes are also of importance, and they will be included in our next set of studies.

Change-point analysis is the main statistical method we employed in this study. The change-point methodology has long been recognized as an important tool for identifying nonstationarities in various scientific phenomena, including climatic and environmental factors. The detection of change points in time series of physical and scientific phenomena is a prominent area of interest in data analysis. The statistical methodology, often addressed as change-point methodology, aims to uncover both the number and positions of temporal transitions in the statistical properties of data. The time series data analyzed in such studies can come in a variety of forms, including univariate or multivariate, independent or dependent, discrete or continuous, low-dimensional or high-dimensional, etc. For a comprehensive understanding of the methodology, one may look at, among others, the monograph by Csörgő and Horváth [17] and also the review article of Jandhyala et al. [18]. Detecting multiple change points in time series is of great interest in most applications. Numerous methodologies have been proposed in the recent change-point literature, including binary segmentation (BS) [19], optimal partitioning (OP) [20], pruned exact linear time (PELT) [21,22], the Bayesian information criterion (BIC) [23], Bayesian approach [24], the minimum description length (MDL) [25,26], graph-based approach [27], and those based on functional time series analysis [28,29]. Recent literature has contributed to advancing the field further. Cheung et al. [30] and Aue and Van Delft [31] have explored recent developments and applications in multiple change-point detection for time series of networks. Additionally, Kaul et al. [32], Gösmann et al. [33], and others have provided insights into detecting changes in high-dimensional time series. Most recently, Fotopoulos et al. [34] and Paparas et al. [35] described methods of both the detection and estimation of change points in gamma- and Poisson-distributed data, respectively. While the basis for their likelihood ratio-based detection method may be found in the work of Csörgő and Horváth [17], the basis for maximum likelihood-based change-point estimation may be found in the works of Jandhyala and Fotopoulos [36] and Fotopoulos et al. [37].

In this article, we first implement the BIC [23] method to identify the number and locations of change points in any given dataset. Since the method is heuristic, the change points identified using this method lack a statistical significance associated with them. For this reason, we implemented the likelihood ratio test to test the statistical significance of each change point identified via the BIC. For the change points that were found to be statistically significant, we constructed confidence intervals, implementing the method of maximum likelihood estimation (mle) [34,35]. A subsequent implementation of the PELT method [21,22] yielded mostly the same statistically significant change points identified through the BIC, and thus, we are confident that the change points identified in this article are quite robust.

Lastly, one may wish to know whether it is appropriate to identify changes in climatic variables through abrupt change-point modeling. For this, one merely needs to look at the vast majority of change-point models that have been applied in the area of meteorology and climatology for the past three to four decades through review articles such as those of Beaulieu et al. [38], Jandhyala et al. [18], Reeves et al. [39], and, most recently, Lund et al. [40]. Also, Pitman and Stouffer [41] focus fully on discussing the relevance of abrupt change models for climate and climate modeling. The recent expository article of Lund and Shi [42] is also an excellent source for understanding the importance of change-point models for identifying changes in climatic variables.

The organization of this paper is as follows. We begin by presenting the following in detail: the study area in Section 2.1, the data sources in Section 2.2, and the statistical methods in Section 2.3. The results from statistical analysis and a detailed discussion of the results, including the change points identified, are presented in Section 3. Section 4 ends the paper with some conclusions.

2. Data and Methods

2.1. Study Area

The Great Lakes region constitutes a multifaceted and dynamic system exerting a considerable influence over North America’s climate and environment. Its susceptibility to diverse atmospheric and oceanic phenomena significantly shapes temperature, precipitation, and snowfall patterns. Understanding the variability and changes in these climatic variables is essential to assessing the repercussions and risks of climate change on the hydrological cycles, ecosystems, and socio-economic sectors within the region. This study delved into annual temperature, precipitation, and snowfall records obtained from six stations located on the rims of Lake Erie and Lake Michigan. These climatic variables have relationships with large-scale climate patterns represented by the North Atlantic Oscillation (NAO) and the Pacific Decadal Oscillation (PDO) indices, which are known to influence the climate variability in the region. The primary objectives of this study were to identify potential change points in climatic variables and briefly explore the possible reasons for these changes.

As described in the introduction, six ASOS stations were selected for inclusion in this study—Buffalo Niagara International Airport (BUF), Erie International Airport (ERI), and Cleveland Hopkins International Airport (CLE), which are located on the southern rim of Lake Erie, and Chicago Midway Airport (MDW), Milwaukee Mitchell Airport (MKE), and Green Bay A S International Airport (GRB), which are located on the western rim of Lake Michigan. Figure 1 displays these six stations on a geographical map. For the purposes of this article, hereafter, we refer to these stations as Buffalo, Erie, Cleveland, Chicago, Milwaukee, and Green Bay.

2.2. Data Sources

For each station, the Global Summary of the Year (GSOY) data files provided by the National Center for Environmental Information (NCEI) of the National Oceanic and Atmospheric Association (NOAA) contain quality-controlled annual summaries of several climatic and meteorological variables computed from stations in the Global Historical Climatology Network. From among all the available variables at each station, we chose nine specific climatic variables representing temperature, precipitation, and snowfall for this study: PRCP1—the number of days in a year with ≥1 inch of precipitation; SNOW1—the number of days in a year with ≥1 inch of snowfall; TMAX32—the number of days in a year with a maximum temperature ≤32 °F; and TMAX90—the number of days in a year with a maximum temperature ≥90 °F. We also chose five continuous variables: PRCP—the total annual precipitation (mm); SNOW—the total annual snowfall (mm); TAVG—the average annual temperature (°C); TMAX—the average annual maximum temperature (°C); and TMIN—the average annual minimum temperature (°C). Data on annual averages for each factor were computed from equally weighted monthly data, with no assigned weights based on the number of days in each month. Annual values were designated as missing if one or more data months during a year were missing. The data spanned through the years of 1939 to 2023, except for Chicago, for which the data began from the year 1942. Generally, data that are released by the NOAA are quite clean, and we did not have to take any data-preparatory steps in proceeding with the data analysis of any of the variables. The only concern was about missing observations. Even in this regard, the missing data were minimal, ranging from a minimum of 0% to a maximum of 1.5% of the number of data points for any given variable.

Table 1. List of the six stations, which are located on the rims of Lake Erie and Lake Michigan.

Location	Geographic Region	Years of Data	% Missing
Buffalo	Lake Erie	1939–2022	0.27
Erie	Lake Erie	1939–2022	1
Cleveland	Lake Erie	1939–2022	0
Chicago	Lake Michigan	1942–2022	1.5
Milwaukee	Lake Michigan	1939–2022	0.66
Green Bay	Lake Michigan	1892–2022	0

depicts the six stations, the region they are located in, the years for which data are included in this study, and the % of data missing at each station. Whenever a data point was missing, we estimated it with the average of the previous five data points.

2.3. Statistical Methods

In this study, we employed change-point analysis methods to identify change points in each of the nine climatic variables. The methodology was developed to detect and identify unknown locations of one or more time points, known as change points, at which given data over time might have changed in their statistical properties. These alterations could include sudden shifts in mean or variance, or changes in the underlying distribution or model of the data. The goal of change-point analysis is to estimate the number of change points and their locations, along with characterizing the data segments before and after each change point. This methodology has found widespread applications in various fields, such as quality control, signal processing, bioinformatics, and economics. In climatology, change-point analysis aids in identifying and understanding the causes and effects of climate variability and change, such as natural cycles, anthropogenic influences, extreme events, and regime shifts [38,43,44].

We applied two methods to detect and identify multiple change points in the climatic variables, namely the BIC [23] and PELT [21,22]. We considered these two methods mainly to obtain greater confidence about the change points that were identified. It turned out that most of the change points identified were common to both the BIC and PELT, and only a few belonged to just one of the methods. Hence, we decided to present results on change points that were identified by implementing the BIC method only.

Both methods are based on the idea of minimizing a penalty cost function that maintains a balance between the fit of the data and the number of change points. The methods assume that data can be modeled using a parametric distribution with a set of unknown parameters that may change at unknown time points. The goal is to estimate both the number and the locations of the change points by minimizing a cost function. For both methods, the cost function is the negative log-likelihood of the data segment, which is proportional to the sum of squared residuals (SSR). The BIC method of Bai and Perron [29] guarantees an exact solution with a linear computational cost. Here, the penalty term is chosen based on the Bayesian information criterion (BIC). The PELT method uses a dynamic programming algorithm to find the optimal solution to the minimization problem, in which the penalty term is chosen according to the modified Bayesian information criterion (MBIC) and is derived as $\beta ={\displaystyle \frac{\mathit{log}(n)}{2}}$, where n is the number of data points. The PELT method also guarantees a computational cost that is linear in the number of data points.

Based on goodness of fit testing, we found the continuous variables to follow either the normal or the gamma distribution, whereas the discrete variables followed the Poisson distribution. We applied both the BIC and PELT methods in all three of these cases. In doing so, we imposed a minimum gap of 20 years between any two consecutive change points, mainly to avoid identifying spurious and excessive changes.

The BIC and PELT methods do not provide statistical significance (p-value) associated with any of the identified change points. In such a case, one wonders how many of the change points identified using either method are truly statistically significant. For this purpose, we performed the likelihood ratio change detection test [17] to assess the statistical significance of each of the detected change points. For the purposes of completeness, we briefly present below the likelihood ratio change detection statistic.

In order to introduce the likelihood ratio statistic, let $Y_1,Y_2,\dots ,Y_n$ be a time series of independently distributed random variables, each with a probability density function (pdf) (or probability mass function, pmf, in the discrete case)$f\left(y;\mathit{\theta }\right)$, where $\mathit{\theta }=\left({\theta }_1,{\theta }_2,\dots ,{\theta }_p\right)$ is a p-dimensional parameter, such that $\theta \in \mathbf{\Theta }$. Suppose a single change-point is detected within the data $Y_{j_1+1},Y_2,\dots ,Y_{j_2}$ using either of the two methods, BIC or PELT. Then, to test for the true presence of this change, we begin by formulating the null hypothesis that there is no change point within $Y_{j_1+1},Y_2,\dots ,Y_{j_2}$ and the alternative hypothesis that there is a change point at an unknown time within $Y_{j_1+1},Y_2,\dots ,Y_{j_2}$. The test is based on the ratio of the likelihoods of the data under the null and alternative hypotheses, and it is given as follows:

$$\Lambda ={\mathit{max}}_{j_1+1\le t\le j_2-1}{\displaystyle \frac{{\mathit{sup}}_{\mathit{\theta }\in \Theta }{\prod }_{i=j_1+1}^{t}f\left(Y_i,\mathit{\theta }\right){\mathit{sup}}_{\mathit{\theta }\in \Theta }{\prod }_{i=t+1}^{j_2}f\left(Y_i,\mathit{\theta }\right)}{{\mathit{sup}}_{\mathit{\theta }\in \Theta }{\prod }_{i=j_1+1}^{j_2}f\left(Y_i,\mathit{\theta }\right)}}$$

Under the null hypothesis, the test statistic has an asymptotic distribution that depends on the distribution of the data. We used this result to compute the p-values for the change points identified via the BIC or PELT, and we considered a point to be statistically significant if the p-value was less than 0.05. For each of the change points found to be statistically significant, it was important to find the confidence interval estimate associated with that unknown point of time.

Each change point found to be statistically significant via the likelihood ratio test may also be viewed as the maximum likelihood estimate (mle) of the unknown change point. Thus, we could apply the asymptotic distribution of the change-point mle as derived by Jandhyala and Fotopoulos [36] in order to compute the confidence interval estimate of any size. Here, we shall first present the asymptotic distribution of the change-point mle for the general case. The computational methods for finding the asymptotic distribution of the change-point mle for the special cases of normal, gamma, and Poisson distributions are presented in Appendix A. As per the change-point set-up, there exists an unknown change-point, $\tau \in \left\{\mathrm{1,2},\dots ,n-1\right\}$, such that

$$Y=\left\{\begin{array}{c}{Y}_t~f\left(y;{\mathit{\theta }}_0\right),t\in \left\{1,\dots ,\tau \right\}\\ Y_t~f\left(y;{\mathit{\theta }}_1\right),t\in \left\{\tau +1,\dots ,n\right\},\end{array}\right.$$

${\mathit{\theta }}_0=\left({\theta }_{10},{\theta }_{20},\dots ,{\theta }_{p0}\right)$is the parameter before a change, ${\mathit{\theta }}_1=\left({\theta }_{11},{\theta }_{21},\dots ,{\theta }_{p1}\right)$ is the parameter after the change, and ${\mathit{\theta }}_0\ne {\mathit{\theta }}_1$. Let ${\widehat{\tau }}_n$ be the maximum likelihood estimator (mle) of $\tau$. Then, as shown by Fotopoulos et al. [37], the centralized change-point mle ${\xi }_n:={\widehat{\tau }}_n-{\tau }_n$ can be expressed as follows:

$${\xi }_n={\mathit{argmax}}_{t\in \left\{-{\tau }_n+1,\dots ,n-{\tau }_n+1\right\}}{C}_n\left(t\right)$$

$$C_n\left(t\right){=}_D\left\{\begin{array}{cc}{S}_{1:\left|t\right|}={\displaystyle \sum _{k=1}^{\left|t\right|}}{X}_k={\displaystyle \sum _{k=1}^{\left|t\right|}}\mathit{ln}{\displaystyle \frac{f\left(Y;{\mathit{\theta }}_1\right)}{f\left(Y;{\mathit{\theta }}_0\right)}},\hfill & t=-{\tau }_n+1,\dots ,-1\\ 0,\hfill & t=0\\ {\displaystyle S_{1:t}^{*}}={\displaystyle \sum _{k=1}^t}{\displaystyle X_k^{*}}={\displaystyle \sum _{k=1}^{\left|t\right|}}\mathit{ln}{\displaystyle \frac{f\left(Y;{\mathit{\theta }}_0\right)}{f\left(Y;{\mathit{\theta }}_1\right)}},\hfill & t=1,\dots ,n-{\tau }_n.\end{array}\right.$$

In the above, the random walks S and S^* are independent of each other. Per Fotopoulos et al. [37], the asymptotic form of ${\xi }_n$, denoted as ${\xi }_{\infty }$, has the following distributional form:

$$P({\xi }_{\infty }=j)\cong \left\{\begin{array}{cc}{e}^{-B}\left[q_{-j}-{\int }_{0+}^{\infty }\left\{1-G^{*}(x)\right\}d{u}_{-j}(x)\right],\hfill & j<0\\ e^{-B-B^{*}},\hfill & j=0\\ e^{-B^{*}}\left[{q^{*}}_j-{\int }_{0+}^{\infty }\left\{1-G(x)\right\}d{u^{*}}_j(x)\right],\hfill & j>0\end{array}\right.$$

Let us define the distribution function of the overall maxima as $G\left(X\right)=P\left(M\le x\right)$and $G^{*}\left(X\right)=P\left(M^{*}\le x\right)$ for $x\ge 0$. We let $M_j=\underset{0\le i\le j}{max}{S}_i$ and ${M_j}^{*}=\underset{0\le i\le j}{max}{S_i}^{*}$be the maxima of the first n partial sums, and we let $M=\underset{n\ge 0}{max}{S}_n$ and $M^{*}=\underset{n\ge 0}{max}{S_n}^{*}$be the total maxima. Also, we let $q_n=P\left({T_X}^{\le }>n\right)$ and $u_n\left(x\right)=P\left(T_X^{\le }>n,S_n\in 0,x\right)$ for $n\ge 0, x\ge 0$ and $B=\sum b_n/n$, where $b_n=P(S_n>0)$, for $n\ge 1$ and $B^{*}=\sum {b_n}^{*}/n$, where ${b_n}^{*}=P({S_n}^{*}>0)$, for $n\ge 1$. Based on the computational methods elaborated in the Appendix A, the asymptotic distribution for each change-point mle was computed depending upon the underlying distribution being normal, gamma, or Poisson.

3. Results from Statistical Analysis and Discussion

The computational methods elaborated in Section 2.3 above enabled us to first implement both BIC and PELT methods to identify change points in each of the data series considered in the article. Subsequently, the likelihood ratio statistic was applied to test for the statistical significance of each change point identified via the BIC and PELT methods. However, since both methods yielded almost the same change points, we present the results of applying only the BIC method. Finally, asymptotic distributions were computed for each statistically significant change point from the BIC, and these enabled us to compute the associated 95% confidence interval estimates. As a follow-up, we also provide estimates of the amount of change in the mean whenever a significant change was identified. However, when a change point identified via the BIC was not statistically significant, we left such point estimates alone without computing the 95% confidence interval or the amount of change. The computed point estimates from the BIC, together with the p-values from the likelihood ratio test, the 95% confidence intervals in the form [point estimate $\pm$ number of years] only when p-value < 0.05, and the amount of change are presented in the following:

Table 2. Change points in mean according to the BIC method for the variables PREC1, SNOW1, TMAX32, and TMAX90: change-point estimate ± width of approximate 95% confidence interval when the change detection test p-value < 0.05 (*), followed by the [amount of change in the mean] or only the point estimate when p-value > 0.05.

	Climatic Variable	PRCP1	p-Value	SNOW1	p-Value	TMAX32	p-Value	TMAX90	p-Value
City
Buffalo		1975 ± 9 [2.36]	0.02 *	1979	0.56	1985 ± 7 [−7.96]	0.01 *	1964 ± 4 [−3.42]	0.00 *
Erie		1969	0.76	1966 ± 8 [5.46]	0.01 *	1985 ± 6 [−8.98]	0.00 *	1965 ± 4 [−2.90] 1990 ± 4 [2.70]	0.00 * 0.00 *
Cleveland		1982 ± 9 [2.30]	0.04 *	1971	0.26	1985 ± 8 [−6.96]	0.01 *	1964 ± 2 [−8.34]	0.00 *
Chicago		1981	0.16	1966 1991	0.27 0.78	1985 ± 7 [−7.49]	0.01 *	1966 1994	0.25 0.79
Milwaukee		1975	0.07	1969	0.40	1985 ± 4 [−11.93]	0.00 *	1964 1991	0.32 0.91
Green Bay		1977	0.18	1970 1996	0.11 0.37	1986 ± 6 [−11.21]	0.00 *	1966 1991	0.65 0.20

for the variables PREC1, SNOW1, TMAX32, and TMAX90;

Table 3. Change points in the mean according to the BIC method for the variables PRCP and SNOW: change-point estimate ± the width of the approximate 95% confidence interval when the change detection test p-value < 0.05 (*), followed by the [amount of change in the mean] or only the estimate when the change detection test p-value > 0.05.

	Climatic Variable	PRCP	p-Value	SNOW	p-Value
City
Buffalo		1974 ± 10 [128.61]	0.02 *	1969	0.36
Erie		1971 ± 11 [134.45]	0.00 *	1963 ± 7 [773.31]	0.02 *
Cleveland		1971 ± 9 [141.76]	0.00 *	1971 ± 15 [317.38]	0.00 *
Chicago		1971 ± 7 [142.98]	0.00 *	1966 1991	0.81 0.90
Milwaukee		1971 ± 8 [153.81]	0.00 *	1970 ± 12 [221.21] 1995 ± 8 [−103.04]	0.00 * 0.00 *
Green Bay		1991 ± 11 [96.59]	0.00 *	1970 ± 3 [310.59] 1995 ± 14 [76.51]	0.00 * 0.00 *

for the variables PRCP and SNOW; and

Table 4. Change points in the mean according to the BIC method for the variables TAVG, TMAX, and TMIN: change-point estimate ± the width of the approximate 95% confidence interval when the change detection test p-value < 0.05 (*), followed by the [amount of change in the mean] or only the estimate when the change detection test p-value > 0.05.

	Climatic Variable	TAVG	p-Value	TMAX	p-Value	TMIN	p-Value
City
Buffalo		1989 ± 8 [0.74]	0.02 *	1986	0.20	1989 ± 6 [0.96]	0.00 *
Erie		1985 ± 6 [1.06]	0.00 *	1985 ± 7 [0.94]	0.01 *	1985 ± 5 [1.18]	0.00 *
Cleveland		1989 ± 8 [0.77]	0.04 *	1989	0.40	1989 ± 5 [1.11]	0.00 *
Chicago		1985 ± 7 [0.96]	0.01 *	1985	0.20	1985 ± 4 [1.33]	0.00 *
Milwaukee		1986 ± 3 [1.32]	0.00 *	1986 ± 7 [0.98]	0.01 *	1985 ± 2 [1.66]	0.00 *
Green Bay		1986 ± 9 [0.88]	0.02 *	1986 ± 9 [0.84]	0.01 *	1989 ± 10 [0.93]	0.04 *

for the variables TAVG, TMAX, and TMIN, respectively.

Also, for purposes of visually understanding the changes, statistically significant change points, along with the corresponding data, are plotted for each data series and for each city. For the purposes of completeness, we also include plots of the rest of the climatic variables for which no statistically significant change point was found. In such cases, we fit just the common mean and denote it as a dashed line. The respective plots are presented in the following: Figure 2 for Buffalo, Figure 3 for Erie, Figure 4 for Cleveland, Figure 5 for Chicago, Figure 6 for Milwaukee, and Figure 7 for Green Bay.

We shall now begin our discussion of the results. For this purpose, we find it convenient to summarize the information in Table 2, Table 3 and Table 4 more concisely, and thus, we present in

Table 5. Summary of changes in the precipitation variables, namely PRCP1—the number of days in a year with ≥1 inch of precipitation, and PRCP—the total annual precipitation (mm). Here, * implies that a change in the corresponding year is statistically significant, (#) implies the corresponding variable is a counting variable, $\uparrow$ implies an increase in a change.

	Lake City	Lake Erie			Lake Michigan
Precip		Buffalo	Erie	Cleveland	Chicago	Milwaukee	Green Bay
PRCP1 (#)		1975 *$\uparrow$	1969	1982 *$\uparrow$	1981	1975	1977
PRCP (mm)		1974 *$\uparrow$	1971 *$\uparrow$	1971 *$\uparrow$	1971 *$\uparrow$	1971 *$\uparrow$	1991 *$\uparrow$

a summary of changes in precipitation, in

Table 6. Summary of changes in the two snowfall variables, namely SNOW1—the number of days in a year with ≥1 inch of snowfall, and SNOW—the total annual snowfall (mm). Here, * implies that a change in the corresponding year is statistically significant, (#) implies the corresponding variable is a counting variable, $\uparrow$ implies an increase in a change, and $\downarrow$ implies a decrease in a change.

	Lake City	Lake Erie			Lake Michigan
Snow		Buffalo	Erie	Cleveland	Chicago	Milwaukee	Green Bay
SNOW1 (#)		1979	1966 *$\uparrow$	1971	1966 1991	1969	1970 1996
SNOW (mm)		1969	1963 *$\uparrow$	1971 *$\uparrow$	1966 1991	1970 *$\uparrow$ 1995 *$\downarrow$	1970 *$\uparrow$ 1995 *$\uparrow$

a summary of changes in snowfall, and in

Table 7. Summary of changes in the five temperature variables, namely the following: TMAX32—the number of days in a year with a maximum temperature ≤32 °F; TMAX90—the number of days in a year with a maximum temperature ≥90 °F; TAVG—the average annual temperature (°C); TMAX—the average annual maximum temperature (°C); and TMIN—the average annual minimum temperature (°C). Here, * implies that a change in the corresponding year is statistically significant, (#) implies the corresponding variable is a counting variable, $\uparrow$ implies an increase in a change, and $\downarrow$ implies a decrease in a change.

	Lake City	Lake Erie			Lake Michigan
Temp.		Buffalo	Erie	Cleveland	Chicago	Milwaukee	Green Bay
TMAX32 (#)		1985 *$\downarrow$	1985 *$\downarrow$	1985 *$\downarrow$	1985 *$\downarrow$	1985 *$\downarrow$	1986 *$\downarrow$
TMAX90 (#)		1964 *$\downarrow$	1965 *$\downarrow$ 1990 *$\uparrow$	1964 *$\downarrow$	1966 1994	1964 1991	1966 1991
TAVG (°C)		1989 *$\uparrow$	1985 *$\uparrow$	1989 *$\uparrow$	1985 *$\uparrow$	1986 *$\uparrow$	1986 *$\uparrow$
TMAX (°C)		1986	1985 *$\uparrow$	1989	1985	1986 *$\uparrow$	1986 *$\uparrow$
TMIN (°C)		1989 *$\uparrow$	1985 *$\uparrow$	1989 *$\uparrow$	1985 *$\uparrow$	1985 *$\uparrow$	1989 *$\uparrow$

a summary of changes in temperature.

Before we delve into details of the nature of changes, it is clear beyond doubt from Table 5, Table 6 and Table 7 that significant changes occurred in precipitation, snowfall, and temperature at all six of the rim stations, and it is important for meteorologists to make note of this significant observation. Truly, it is meteorologists who can access the depths of scientific phenomena and discover the reasons that have contributed to the significant changes identified in this study.

From Table 5, it is clear that the total annual precipitation (PRCP) increased consistently around the year 1971 at all the stations except for Green Bay, where the increase occurred in 1991. However, significant increases in the number of days with precipitation (PRCP1) occurred only in Buffalo in 1975 and in Cleveland in 1982. It is important to note that all significant changes in both PRCP and PRCP1 were on the increasing side for precipitation. Such consistently increasing changes did not occur for either snowfall or temperature. In this sense, changes in precipitation stand out and require further attention.

With respect to changes in snowfall, Table 6 shows that snowfall at the rim stations of Lake Erie increased consistently, whereas changes in snowfall at rim stations of Lake Michigan did not show similarly consistent increases. The presence of two changes in snowfall (SNOW) occurred at Green Bay with increases in both 1970 and in 1995, whereas in Milwaukee, there was initially an increase in 1970, and then the annual snowfall decreased significantly in 1995. With regard to changes in the number of days in a year with ≥1 inch of snowfall (SNOW1), significant increases occurred only at Erie in 1966. Otherwise, no significant changes occurred in SNOW1 at other stations of Buffalo, Cleveland, Chicago, Milwaukee, or Green Bay.

Table 7 shows that significant increases in the continuous temperature variables—the average annual temperature (TAVG), average annual maximum temperature (TMAX), and average annual minimum temperature (TMIN)—occurred at all six stations. It is noteworthy that most of these increases occurred in or around the year 1985. The only exceptions were the stations of Buffalo and Cleveland, where significant increases in TAVG and TMIN occurred around 1989, not too long after 1985. In contrast, changes in the number of days with a maximum temperature $\le$32 °F (TMAX32) decreased consistently at all six stations, also around the same year of 1985. In terms of changes in the number of days with a maximum temperature $\ge$90 °F (TMAX90), all significant changes occurred only at the three rim stations of Lake Erie, and no significant changes were found at the rim stations of Lake Michigan. There were two significant changes at Erie Station, with the first change being a decrease in TMAX90 around 1965, while the second change led to an increase in TMAX90 around 1990.

As for the extent of changes in precipitation, snowfall, and temperatures, we limit our discussion to the continuous variables only. From Table 3 and Table 4, we note that the highest increase in precipitation (PRCP) can be seen to have occurred in Milwaukee (153.81 mm), followed by Chicago (142.98) and Cleveland (141.76). Since Cleveland, Chicago, and Milwaukee are geographical neighbors to one another, anticipating a similar trend, one would have thought Green Bay would have undergone an even higher increase in precipitation than Milwaukee. However, it should be noted that the higher increases that occurred in Cleveland, Chicago, and Milwaukee all happened in 1971, whereas the lowest increase in precipitation that occurred at Green Bay happened much later, in 1991. The greatest extent of change in snowfall was an increase to the tune of 773.31 mm, which happened in 1963 at Erie. Next, Cleveland exhibited a very large increase (317.38 mm) in snowfall in 1971, whereas Chicago underwent no significant increase in snowfall. Milwaukee exhibited first an increase in 1970 (221.21 mm) and then a decrease in 1991 (−103.04 mm), while the lowest increase (76.51 mm) in snowfall occurred in the year 1995 at Green Bay.

Let us move on to the temperature factors of TAVG, TMAX, and TMIN. The highest increase in temperature (1.66 °C) in the TMIN factor occurred at Milwaukee in 1985, and this was followed by the next highest increase of 1.33 °C in the same variable, TMIN, which occurred in the neighboring Chicago in the same year of 1985. The third highest increase of 1.18 °C in TMIN occurred again in 1985 in the city of Erie. It should be mentioned that, with an increase of 1.11 °C in TMIN, Cleveland was not far behind. A somewhat similar pattern can be seen for the extent of changes in TAVG. The highest increase of 1.32 °C happened at Milwaukee, followed by an increase of 1.06 °C in the city of Erie, and the increase at Cleveland stood at 0.77 °C. As for TMAX, the increases were significant only at Milwaukee (0.98 °C), Erie (0.94 °C), and Green Bay (0.84 °C), and these changes happened essentially in the year 1986.

It is also important to discuss the 95% confidence intervals that are presented in Table 2, Table 3 and Table 4. Once again, we limit our discussion here to the continuous variables only, and hence, we limit ourselves to Table 3 and Table 4 only. Since all the confidence intervals were computed at the same 95% level, the corresponding widths were all directly comparable. From the viewpoint of the shortest 95% confidence intervals, we note that TMIN and TAVG at Milwaukee, with intervals of $\pm$2 and $\pm$3, respectively, were the shortest of all the confidence intervals. The next shortest confidence interval of $\pm$4 occurred for TMIN at Chicago. The next shortest interval of $\pm$5 also occurred for TMIN at both Cleveland and Erie. It is notable that almost all the confidence intervals for PRCP and SNOW were uniformly high, with the largest intervals occurring at Cleveland ($\pm$15, SNOW) and Green Bay ($\pm$14, SNOW). From these confidence intervals, it may be concluded that change-point estimates for PRCP and SNOW are less reliable compared to those for the temperature variables.

One may notice that there were consistent decreases in 1985 at all six stations in the number of days with a maximum temperature $\le$32 °F (TMAX32) and also uniform increases that occurred at almost all stations in TAVG, TMAX, and TMIN in and around the same year of 1985. The observed decreases in the number of days not exceeding 32 °F are consistent with increases in the continuous variables TMIN, TAVG, and TMAX since a general shift in the distributions towards warmer temperatures would naturally decrease the frequency of especially cold temperatures. However, the significant decrease in TMAX90 in the rim stations of Lake Erie in and around 1965, and then a significant increase at Erie Station in and around 1990, requires a careful understanding and explanation. Another observation of interest is that there were fewer changes in the four discrete variables than there were in the five continuous variables. It could be argued that the information content in the discrete variables was much less compared to the information contained in the continuous variables, and this may have led to greater sensitivity in the change detection statistic. In this sense, the fewer changes to the discrete variables could be viewed as being conservative.

Instead of the change-point approach adopted in this article, climatologists also study changes in climatic variables alternatively through simple linear trend models (e.g., see Lai and Dzombak [45] and Isaac and Wijngaarden [46]). It then becomes relevant to know which of the two models fits better for the climatic data considered in this study. The Akaike Information Creterion (AIC) is a widely applied method for such model selection problems. While there are other information criteria-based methods, such as the Bayesian information criterion (BIC) and so on, the AIC method (see Akaike [47]) is one of the most frequently adopted methods for model-selection problems. For any given model, the AIC is given as follows:

AIC = −2lnL + 2k

where L is the likelihood function, and k is the number of parameters in the model under consideration. Upon computing the AIC number for all models, one selects the model with the smallest AIC number as the best model.

In our case, we have two competing models, the change-point model and the simple linear trend model. We computed the AIC numbers for both models for each of the 54 combinations of stations and climatic variables, and the resulting AIC numbers are presented in

Table 8. Akaike Information Criterion (AIC) number for both the change-point model and the simple linear trend model for all nine climatic variables at all six stations.

AIC Numbers for Change-Point Model
Station	PRCP1	SNOW1	TMAX32	TMAX90	PRCP	SNOW	TAVG	TMAX	TMIN
Buffalo	375.01	549.94	663.10	481.09	1060.49	1293.99	190.60	200.87	195.07
Cleveland	398.00	501.24	656.72	606.06	1098.47	1213.51	205.32	221.61	205.12
Erie	433.63	584.48	669.48	454.37	1104.22	1351.99	212.88	217.56	223.54
Milwaukee	400.17	510.57	671.08	561.16	1077.86	1261.25	209.04	218.47	213.87
Chicago	406.38	460.89	634.63	621.49	1063.15	1191.81	195.04	206.16	199.98
Green Bay	384.13	508.59	660.08	537.67	1046.48	1240.95	218.38	213.31	245.36
AIC Numbers for Simple Linear Model
Buffalo	386.16	551.85	669.53	488.71	1065.51	1298.75	194.19	204.92	193.64
Cleveland	399.03	503.08	659.92	616.55	1096.41	1223.94	212.85	224.77	213.38
Erie	436.15	584.38	675.10	460.18	1107.72	1352.96	225.44	224.75	240.65
Milwaukee	402.48	509.29	678.79	560.06	1080.82	1271.57	223.06	228.38	228.30
Chicago	406.37	462.50	638.32	620.86	1058.65	1197.62	201.84	208.94	207.53
Green Bay	386.64	509.58	669.73	538.58	1038.73	1253.98	226.73	223.00	251.94

below.

The minimum AIC number criterion clearly shows that except for 8 situations, the change-point model is the preferred model for all the remaining 46 combinations of stations and climatic variables. Hence, the AIC criterion prefers the change-point model in an overwhelming number of cases thus showing the relevance of change-point approach for modeling changes in the climatic variables at the 6 lake rim stations.

Finally, It is important to check for the validity of the assumption of independence of the time series over time. At first glance, it may appear that there is the presence of autocorrelations in the original data series. However, once we fit the change-point model and thus account for the change points, then the corresponding residuals do not show the presence of autocorrelations any further. We have verified this through ACF (auto-correlation function) and PACF (partial autocorrelation function) plots of both the original data, and also of the residuals derived from the change-point model. To maintain brevity, we do not include the ACF and PACF plots in the article.

4. Conclusions

This study found significant changes in precipitation, snowfall, and temperatures at all six stations located on the rims of Lake Erie and Lake Michigan. Before concluding this article, it is important to compare the results of this study to the results of similar studies in the literature. However, while searching, it became difficult to find studies that were quite similar in scope. Nevertheless, there have been some studies that have assessed changes in the overall Great Lakes area and also in the midwestern US, and we compared our results with those of such studies. For the purposes of comparison, we utilized the amounts of changes at all six stations in Table 2, Table 3 and Table 4 and computed the following: the average amount of change for TAVG = 0.96 °C, the average percentage of change for PRCP = 16%, and the average percentage of change for SNOW = 17%. Wuebbles et al. [48] prepared an assessment of the impact of climate change on the Great Lakes, and in that assessment, they reported that, from 1901–1960 to 1985–2016, the annual mean temperature in the Great Lakes basin increased by about 1.6 °F or 0.89 °C. While their reported change refers to the entire Great Lakes basin, in our study of six stations on the rims of Lake Erie and Lake Michigan, we found that the annual mean temperature (TAVG) increased by about 0.96 °C. As for precipitation, Wuebbles et al. [48] reported an increase of 10% for the Great Lakes basin, whereas we found that precipitation increased, on average, by about 16% for the two lakes. In the same report, Weubbles et al. [48] reported that the annual snowfall in the Great Lakes basin decreased by 2.25%. In contrast, in this article, we found that the annual snowfall (SNOW) for the two-lake region increased by as much as 17%. It remains to be investigated whether this discrepancy stems from the fact that the change in snowfall in our article pertains to only the two lakes, whereas Weubbles et al. [48] computed the change in annual snowfall for the entire Great Lakes basin. In a subsequent report on the assessment of climate change in Illinois, Weubbles et al. [49] reported that the average daily temperature in Illinois increased by 1–2 °F or 0.56–1.11 °C and that precipitation increased by 5 to 20%. These changes compare well with our estimated changes of 0.96 °C for temperature and 16% for precipitation. In a recent study, Xue et al. [50] reported that the Great Lakes basin is projected to warm by 1.3–2.1 °C, and precipitation by 0–13%, by the mid-21st century. The changes we found in this article do not extend all the way to the mid-21st century, and hence, our changes of 0.96 °C for temperature and 16% for precipitation are not directly comparable.

It is unclear at this point whether the changes identified in this article can be attributed to changes that might have occurred in Lake Erie and Lake Michigan or more generally in the Great Lakes system. Such an association, if not causation, requires an in-depth study of changes in the Great Lakes. A confounding factor that can complicate establishing an association is the fact that the rim stations of Lake Erie are all climatologically downwind, whereas the rim stations of Lake Michigan are upwind. It is important to understand the effect of this climatological difference while establishing an association between climatic changes at rim stations and changes in lakes’ climatic conditions. As pointed out in the introduction, McBean and Motiee [15] analyzed precipitation, temperature, and stream flows of the Great Lakes over the period of 1930–2000, and Briley et al. [16] reported that the Lake Superior surface waters have warmed more rapidly than nearby air temperatures. While these studies are relevant, they are insufficient to make much headway in explaining climatic changes in rim stations based on changes that might be happening in the Great Lakes system. It is also relevant to ascertain whether changes found in the rim stations can also be found in stations that are relatively farther away from the rim of the Great Lakes. In such a scenario, the changes in the rim stations may not necessarily be due to changes in the Great Lakes. Overall, this study raises the need for other, larger follow-up studies, and it is our endeavor to pursue these studies in the future in a timely manner.