This study examines the primary atmospheric controls over winter precipitation variability in the Great Lakes basin and the potential for seasonal prediction. We employ partial least squares (PLS) regression to identify the primary modes of joint variability between winter precipitation over each of the Great Lakes and concurrent anomalies in midlevel atmospheric flow. We find that the first identified pattern (PLS1) is related to El Niño–Southern Oscillation (ENSO), while the other patterns represent unique anomalies in atmospheric flow that govern precipitation gradients over the basin, with limited seasonal predictability. Nonlinearities are found in the relationship between a sea surface temperature (SST)-based index for ENSO and PLS1 with respect to the phase, strength, and type of ENSO event. An examination of the ENSO-related propagating wave train that drives variability of PLS1 precipitation reveals that seasonally lagged tropical Pacific convection, as measured by remotely sensed outgoing longwave radiation (OLR), is more strongly and linearly related to Great Lakes winter precipitation than SST-based ENSO indices. Cross-validated linear regressions based on October OLR signals explain 20%–32% of the out-of-sample precipitation variability in the Great Lakes basin. We conclude with a deeper assessment of the underlying relationship between patterns of OLR anomalies in the western equatorial Pacific and Great Lakes winter precipitation. Results show that precipitation response to El Niño is similar regardless of OLR intensity in the tropical Pacific, but for La Niña events, the precipitation response is stronger under weak tropical OLR anomalies. The potential for further improvements in ENSO-based seasonal forecasts are discussed.
The Laurentian Great Lakes form the largest group of freshwater lakes on Earth and are central to the environment and economy of central North America (Gronewold et al. 2013). The water balance of the Great Lakes is driven by a combination of direct inflow, overlake precipitation, overlake evaporation, and outflow, but regional precipitation and overlake evaporation are ultimately the moisture sources and sinks that drive changes in lake storage (Gronewold and Stow 2014). For extreme high water levels with severe, adverse effects to lakeshore interests, above-average precipitation plays a dominant role (Gronewold et al. 2016). A recent assessment of major coastal flooding events due to high lake levels on Lake Ontario determined that above-average outflows from Lake Erie are important in the development of high water events, and cumulative winter precipitation over the Great Lakes basin sustain these higher rates of upstream discharge over both winter and spring seasons (Carter and Steinschneider 2018). This motivates the present study, which seeks to determine the large-scale climate patterns that drive winter precipitation over the Great Lakes, and the potential for seasonal prediction to inform flood preparedness and water level management of the lake system.
There is a rich literature on the associations between large-scale modes of atmospheric and oceanic variability [e.g., Pacific–North American (PNA), North Atlantic Oscillation (NAO), El Niño–Southern Oscillation (ENSO), Pacific decadal oscillation (PDO), Atlantic multidecadal oscillation (AMO)] and hydrometeorological variables of the Great Lakes, including water levels (Polderman and Pryor 2004; Ghanbari and Bravo 2008; Hanrahan et al. 2010, 2014; Assani et al. 2016), temperature and ice cover (Assel 1998; Assel et al. 2000; Rodionov and Assel 2000; Notaro et al. 2006; Wang et al. 2012; Bai et al. 2012, 2015), cyclone activity (Angel and Isard 1998; Isard et al. 2000; Lofgren and Bieniek 2008), and precipitation (Leathers et al. 1991; Rodionov 1994; Grover and Sousounis 2002; Steinschneider and Lall 2015, 2016a,b). Of the teleconnections studied, wintertime ENSO and its relation to the PNA have received the most attention, because the skill of ENSO forecasts is better than for other teleconnections, and both ENSO and the PNA provide the greatest potential for seasonal prediction of Great Lakes climate. Rodionov (1994) and others (Angel and Isard 1998; Notaro et al. 2006; Hanrahan et al. 2010) showed that during a positive PNA event, flow in the polar jet stream is more meridional, with an enhanced ridge over the Rocky Mountains and trough over the eastern United States, supporting cyclogenesis near Alberta (Alberta clippers) with storm tracks that enter the Great Lakes with limited moisture. Conversely, during a negative PNA event, the jet stream is more zonal, steering storm tracks along the U.S.–Canadian border and inducing cyclogenesis in the “Colorado low” region, which transports warm, moist air from the Gulf of Mexico into the Great Lakes basin. This phenomenon can be predicted 1–6 months in advance based on ENSO and the well-understood teleconnection linking tropical Pacific sea surface temperatures (SSTs) to the midlatitude atmosphere over the Pacific–North American region (Horel and Wallace 1981).
Despite the potential for seasonal prediction, the relationship between SST-based indices of ENSO and winter precipitation averaged across the Great Lakes is noisy (Fig. 1) and thus may be of limited use for water management applications. There are several factors driving variability in this relationship. First, precipitation can vary greatly across the 1200-km west-to-east span of the Great Lakes basin (Clites et al. 2014) and is subject to different aspects of the ENSO teleconnection. In addition, the atmospheric response to ENSO can be nonlinear, obscuring the relationship between ENSO and surface climate. One manifestation of this nonlinearity is the asymmetric patterns of atmospheric flow, temperature, and precipitation that emerge under El Niño and La Niña events, which are related to a zonal shift of strong convection over the equatorial Pacific between the two ENSO phases (Hoerling et al. 1997, 2001; An and Jin 2004).
The atmospheric response can also differ between events within a particular phase of ENSO, leading to variability in observed regional surface climate during different El Niño or La Niña events. Some of this variation is linked to random atmospheric variability (Hoerling and Kumar 1997; Deser et al. 2017) that emerges on intraseasonal time scales in the extratropics and perturbs, degrades, or otherwise alters the teleconnection pattern (Robertson et al. 2015; Muñoz et al. 2015; Moron et al. 2017), inhibiting seasonal prediction. However, there are also structured differences in the spatial pattern, strength, and temporal evolution of ENSO events (Hoerling and Kumar 2002), with unique surface climate responses that may be predictable. For instance, surface climate in certain regions does not scale linearly with the strength of an El Niño event, due to zonal shifts in tropical Pacific convection and the resulting wave train linked to event strength (Hoell et al. 2016). Spatial patterns of SSTs associates with ENSO have also recently been divided into two distinct “flavors,” the central Pacific (CP) and the eastern Pacific (EP) ENSO, which are defined by the longitudinal location of SST anomalies (SSTAs) over the equatorial Pacific (Kao and Yu 2009; Kug et al. 2009; Capotondi et al. 2015; Feng et al. 2017). These two ENSO typologies can excite distinct responses in atmosphere flow and regional climate (e.g., Larkin and Harrison 2005; Kim et al. 2009; Yu et al. 2012), which a simple index for ENSO (e.g., Niño-3.4) cannot capture. Some have argued that any index based on Pacific SSTs may not fully capture the effect of ENSO on midlatitude atmospheric dynamics, which is more directly related to atmospheric heating driven by anomalous deep convection over the tropical Pacific (Chiodi and Harrison 2010). This process is related to but not entirely determined by SSTs. A direct measure of tropical convection, like outgoing longwave radiation (OLR) over the equatorial Pacific, may better capture the tropical–extratropical relationship (Chiodi and Harrison 2013, 2015; L’Heureux et al. 2015).
There are also other modes of oceanic and atmospheric variability that interact with the ENSO phenomenon to influence climate over the Great Lakes basin (Rodionov and Assel 2000). Some of these teleconnections may themselves be predictable. For example, ENSO-related anomalies in winter temperature over the region are modulated by the NAO, leading to nonlinearity in the underlying teleconnections (Notaro et al. 2006; Bai et al. 2012, 2015). While a large fraction of NAO variance is associated with unpredictable atmospheric noise, there is some potential for seasonal predictability based on summer/autumn Atlantic SST anomalies that are forced by NAO-induced turbulent heat exchanges, but then induce a positive feedback to the atmosphere in the following winter (Gastineau and Frankignoul 2015). In addition, lower-frequency modes of climate variability in the Pacific (PDO) and Atlantic (AMO) can influence Great Lakes climate (Ghanbari and Bravo 2008; Hanrahan et al. 2010, 2014; Assani et al. 2016), and can either interfere with the effects of ENSO directly (Fuentes-Franco et al. 2016) or through their effects on other teleconnections (Gastineau et al. 2013).
Based on the literature above, an assessment of predictability of Great Lakes winter precipitation requires that we understand the atmospheric response to interactions between ENSO and other teleconnection patterns, and differences in Great Lakes precipitation response to structured variations in ENSO events, including any intrinsic nonlinearity in those relationships and sensitivities to how ENSO is measured. There has been some effort in this regard for other hydrometeorological variables (e.g., ice cover, temperature) in the Great Lakes basin (Bai et al. 2012, 2015), and precipitation across broad regions of North America (Notaro et al. 2006; Strong and Liptak 2012; Yu and Zou 2013), but a strict focus on patterns of precipitation across each of the Great Lakes and their nuanced relationship to ENSO and other teleconnections is underdeveloped.
This work contributes a detailed investigation of patterns of winter precipitation across the Great Lakes, and their link to prevailing teleconnections in the Pacific and Atlantic basins, with a focus on different characteristics of ENSO. We first employ partial least squares (PLS) regression to identify the primary modes of joint variability between winter precipitation over each of the Great Lakes and concurrent anomalies in midlevel atmospheric flow, and compare the identified patterns to well-known teleconnections. All but the first pattern are only weakly related to predominant atmospheric teleconnections, and instead represent unique anomalies in atmospheric flow that govern precipitation gradients over the basin. The first mode is significantly but nonlinearly related to tropical Pacific SSTs and extratropical height anomalies commonly associated with ENSO events. We further examine how the propagating wave chain associated with ENSO relates to Great Lakes precipitation, and find that seasonally lagged tropical Pacific convection is more directly and linearly related to Great Lakes winter precipitation than other SST-based ENSO indices, although with some important caveats related to the strength of the OLR signal. Other regions of seasonally lagged convection also exhibit strong predictive signals, while SSTs in the extratropics provide little information for seasonal prediction. We use these results to demonstrate the potential to improve seasonal prediction of winter precipitation variability in the Great Lakes basin.
Time series of monthly overbasin precipitation for each of the Great Lakes were collected from the Great Lakes Monthly Hydrologic Database (Hunter et al. 2015), maintained by the NOAA Great Lakes Environmental Research Laboratory (GLERL). These data were cumulated for the winter (DJF) season and then centered and scaled to produce a 64-yr (1951–2014) series of standardized precipitation anomalies for Lakes Superior (SUP), Michigan–Huron (including Georgian Bay; MHG), Erie (ERI), and Ontario (ONT). Michigan–Huron is considered a single lake unit because the two are hydrologically joined through the Straits of Mackinac. To extend the series until 2017, we downloaded gridded monthly precipitation (0.5° × 0.5° resolution) between 1951 and 2017 from the Global Precipitation Climatology Centre (GPCC; Schneider et al. 2015). The GPCC precipitation was also aggregated to winter anomalies and was converted to standardized precipitation anomalies averaged over each lake basin for 2015–17. This approach was taken because the GLERL database provides the highest-quality precipitation data for the Great Lakes, but a strong El Niño event occurred during the winter of 2016.
Monthly means of 500-hPa geopotential heights between 1950 and 2017 were extracted from the National Centers for Environmental Prediction (NCEP)–National Center for Atmospheric Research (NCAR) Reanalysis 1 dataset (2.5° × 2.5° resolution; Kalnay et al. 1996). We also collected NCEP–NCAR Reanalysis 1 average monthly 300-hPa meridional and zonal wind speeds over the same time period. Absolute, DJF-averaged wind speeds were calculated and used to portray the location of the core of the upper-level jet over North America.
Monthly average SSTs (1° × 1° resolution) between 1950 and 2017 were retrieved from the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST1; Rayner et al. 2003), provided by the Met Office Hadley Centre. Monthly OLR data were obtained from the NOAA Interpolated OLR dataset (Liebmann and Smith 1996), with a shorter time span from June 1974 to January 2017 (data missing in 1978) and a spatial resolution of 2.5° × 2.5°. The OLR data act as a proxy for convection, with lower (higher) OLR indicating convection reaching higher (lower) heights.
Several indices of well-known teleconnection patterns linked to North American climate were also collected. A monthly Niño-3.4 index was collected from the Royal Netherlands Meteorological Institute (KNMI), based on the ERSST v5 dataset (Huang et al. 2017). Following Kao and Yu (2009), EP and CP variants of SST-based ENSO indices were calculated using a regression–EOF method based on the HadISST1 data. In this procedure, a regression is calculated between the Niño-4 (Niño-1+2) index and SSTAs for each grid cell between 20°S and 20°N and 120°E and 80°W. The EP (CP) index is then defined as the first principal component of the space–time field of residuals across that tropical Pacific region. Selected atmospheric circulation teleconnections include the PNA, the NAO, the Greenland Block (Hanna et al. 2013), the tropical/Northern Hemisphere (TNH; Mo and Livezey 1986), and both the western Atlantic (WA) and eastern Atlantic (EA; Wallace and Gutzler 1981). Monthly time series of the PNA index, NAO index, and TNH index were downloaded from the NOAA Climate Prediction Center (http://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml). A daily Greenland Blocking Index (GBI) based on the NCEP–NCAR Reanalysis was retrieved from the NOAA Earth System Research Laboratory (ESRL) Physical Sciences Division (https://www.esrl.noaa.gov/psd/data/timeseries/daily/GBI/). The WA and EA indices were calculated using 500-hPa geopotential heights from the NCEP–NCAR Reanalysis 1 database using the approach in Wallace and Gutzler (1981). All indices were seasonally averaged for each winter between 1951 and 2017.
The methods of this work are composed of four primary analyses. First, we use PLS regression to identify the primary modes of joint variability between Great Lakes winter precipitation and midlevel atmospheric pressure. Composites are then used to explore the underlying circulation and oceanic forcing associated with these modes. Based on the results of those analyses, we examine the relationship between ENSO and Great Lakes winter precipitation in more detail, including seasonal forecast skill based on linear and nonlinear regression models using SST and OLR indices. Finally, we explore whether seasonal predictability of Great Lakes winter precipitation can be better understood using a threshold-based approach for characterizing OLR anomalies in the equatorial Pacific. These analyses are described in more detail below.
a. Joint variability between Great Lakes winter precipitation and midatmospheric pressure
We first employ PLS regression to identify the primary modes of joint variability between DJF precipitation over each of the Great Lakes and concurrent 500-hPa height anomalies over the region from 180° to 20°E and from 10° to 80°N. PLS is a statistical approach to model the linear relationships between two sets of variables, by projecting each set in a latent direction along which the projected variables have maximal covariance (Wold 1975). Let () and () be matrices for 500-hPa geopotential heights and precipitation, respectively, where is the number of years, is the number of grid cells of geopotential height, and is the number of lakes. PLS regression finds two optimal weight vectors w and c that maximize the covariance between the projection of onto w and the projection of onto c:
Here, u and v are termed score vectors. The precipitation time series for each individual lake is centered and scaled into a normalized series so that each lake is considered equally in the optimization above, while the columns of have only been centered (i.e., geopotential height anomalies) to emphasize regions of high atmospheric variability.
The optimization problem in Eq. (1) can be solved by either nonlinear iterative partial least squares (NIPALS; Wold 1975) or eigendecomposition (see Höskuldsson 1988). To extract more than one set of loading patterns and associated score vectors, the above procedure needs to be repeated iteratively after removing previously identified patterns from the data. If , are the score vectors of the first PLS component, we regress each column of on and each column of on , respectively, and retain the residuals:
where and are regression coefficients for the first PLS component, and and are residuals independent of the first PLS component. The optimization in Eq. (1) is then repeated using and to extract the second PLS component, and this process is repeated until components are determined.
Each PLS component represents an independent mode of covariability between atmospheric circulation over a region encompassing the jet entry and exit regions of North America and winter precipitation across each of the Great Lakes. Similar to principal component analysis (PCA), the intensity and phase (positive or negative) of the score vectors for geopotential height u and precipitation v (hereafter referred to as PLSg and PLSp, respectively) indicate the strength and direction of a corresponding spatial pattern in those variables, represented by the associated loading vectors w and c. PLS regression provides a compromise between PCA, which identifies linear patterns of maximum variance for one set of variables, and canonical correlation analysis (CCA; see Hotelling 1936), which finds linear patterns of maximum correlation between two sets of variables. By maximizing covariance between sets of variables, PLS regression can better ensure the identification of highly correlated spatial patterns that explain significant variance in each dataset. This feature is particularly pertinent to the (unstandardized) geopotential height used in this study, in order to appropriately weight the atmospheric anomalies in the extratropics that generally have larger variances than in the subtropics and are known to influence jet dynamics that are important for Great Lakes climate (Rodionov 1994; Rodionov and Assel 2000).
b. Diagnostic analysis of modes of joint variability
The PLS method identifies patterns of winter precipitation across the Great Lakes that are linearly related to large-scale patterns of atmospheric anomalies. While useful, these patterns cannot capture nonlinearity in the underlying relationship. To address this shortcoming, we composite North American precipitation, geopotential heights, and upper-level winds for years where the score vectors of each PLS precipitation pattern are anomalously positive (≥1) or negative (≤−1). These composites reveal how the jet location enhances or suppresses precipitation over each of the Great Lakes as part of a larger pattern of winter precipitation across North America. Effectively, we use PLS as a filter to identify patterns of winter precipitation variability across the Great Lakes that are highly related to different modes of atmospheric flow, but then use composites to examine the nonlinear relationships between those precipitation patterns and their atmospheric drivers. For comparison, the PLS patterns and associate composites are compared to those derived using PCA. We also examine the rank-based Spearman correlation between the score vectors of PLS precipitation patterns and ocean–atmospheric teleconnection indices to identify monotonic relationships with well-recognized climate patterns. Finally, we composite September–November (SON) SSTAs that precede both positive and negative phases of each PLS precipitation pattern, to assess whether ocean boundary conditions during the fall provide any potential to seasonally forecast different modes of Great Lakes winter precipitation. The significance of the SSTA composites is assessed through bootstrapping (Efron 1992).
c. Different measures of ENSO and cross-validated seasonal predictability
We will show that the first PLS mode of Great Lakes precipitation projects strongly onto an ENSO signal. To better understand the monthly progression of ENSO-related phenomena that precedes the first PLS mode, we examine rank correlations between DJF PLSp1 scores and monthly SSTs, OLR, and 500-hPa geopotential heights on a gridcell basis from September to February. These maps show the spatiotemporal evolution of ocean temperatures, convection, and the propagating wave chain associated with ENSO as it relates to the first mode of Great Lakes winter precipitation.
Based on these analyses, we identify regions of SSTs and OLR that are most strongly related to Great Lakes winter precipitation at 1–3-month lead times. We then use linear and generalized additive models (GAMs; see Hastie 2017) to develop seasonal predictions of winter precipitation variability in the Great Lakes basin, with a focus on whether forecasts are substantively improved when using OLR over SSTs. GAMs can account for nonlinearity in a relationship by modeling the expected value of PLSp1 as a smooth function of K predictors:
where the functions are set as smoothing splines with order and parameters estimated using generalized cross validation (Wood 2006). A GAM is used only when substantial nonlinearity is observed in the relationship between PLSp1 and a predictor. Otherwise, simple linear relationships are employed. All models are tested in a leave-one-out cross validation procedure, and we report performance statistics between these predictions and each left-out observation.
Based on the intermodel comparison above, we couple the best predictor set with additional SST indices based on the SON composite maps developed in section 3b, and test several seasonal forecast models for each of the Great Lakes in a leave-one-out cross-validation framework. The motivation for this stepwise approach is that the ENSO-based teleconnection, which relates well to the first PLS mode, will apply broadly to all of the Great Lakes and can be optimized (in terms of covariate selection) based on PLS1. The other PLS modes represent more nuanced precipitation gradients across the Great Lakes basin, and their associated lagged SST signals will likely be more applicable to some lakes over others.
d. OLR El Niño and La Niña event-based index
The continuous forecasting model based on OLR in section 3c is similar to the modeling framework presented in L’Heureux et al. (2015). Harrison and Chiodi (2017) recently critiqued this approach, arguing that the underling relationship between tropical Pacific OLR and surface climate across North America and the globe is more binary and based on whether anomalous OLR in specific regions of the equatorial Pacific crosses a threshold that defines an OLR event. The relevant region of OLR depends on the type of ENSO event. For El Niño, Chiodi and Harrison (2013) defined an OLR El Niño event based on peaks within a time series of 30-day moving-average anomalous OLR within the region 5°S–5°N, 160°–110°W. For La Niña, OLR event years are based on peaks in the cumulative number of days from 1 March to 30 November that anomalous daily OLR within the region 5°S–5°N, 150E°–180° exceeds a threshold (260 W m−2) associated with the 90th-percentile clear-sky value (Chiodi and Harrison 2015). Based on these previous approaches, we identify OLR event years for both El Niño and La Niña between 1975 and 2017, as well as non-OLR El Niño and La Niña event years that do not meet these OLR criteria but are considered El Niño/La Niña events based on the NOAA historical definition. Winter (DJF) precipitation, OLR, and 500-hPa geopotential height anomalies are then composited for OLR and non-OLR event years, as is winter precipitation for each of the Great Lakes. Statistical significance of the composite precipitation averages for each of the Great Lakes is assessed for the different non-OLR and OLR event categories using bootstrapping. The goal of this analysis is not to settle the ongoing debate on the underlying form (continuous, binary) of the ENSO teleconnection to global surface climate (L’Heureux et al. 2017; Harrison and Chiodi 2017), but rather to better understand the implications of these two views for seasonal prediction in the Great Lakes basin.
a. PLS patterns, teleconnections, and composites of atmospheric flow
Figure 2 shows the loading vectors [i.e., c in Eq. (1)] and score vectors [v in Eq. (1), also denoted PLSp] for the four PLS precipitation patterns. We do not show the corresponding loading and score vectors for the geopotential height field, as their relationship to the precipitation patterns are analyzed later using composites. The loading vectors in Fig. 2 highlight patterns of winter precipitation across the lakes that relate strongly to patterns of atmospheric flow, and the associated scores reflect the magnitude and direction of these precipitation patterns through time. PLSp1 assigns large, positive weights to all lakes, representing the average precipitation over the entire Great Lakes basin, although Lake Ontario is weighted somewhat less in this spatial average. As will be shown later, this is likely because climate over Lake Ontario is less influenced by ENSO than the other lakes (Bai et al. 2012). PLSp2 represents a dipole of precipitation across Lake Superior and Lakes Erie and Ontario, while the third pattern is strongly determined by the contrast between precipitation over Lakes Michigan–Huron and Erie. The last mode is characterized by a dipole between Lake Erie and Lake Ontario, with little influence from the upper lakes.
Table 1 shows the Spearman correlation coefficient between the PLSp score vectors and a variety of teleconnection indices. Of all of the PLSp patterns, PLSp1 has the strongest monotonic relationships to the teleconnection indices considered here. The correlation coefficients (CCs) between PLSp1 and DJF Niño-3.4 and the PNA are −0.33 and −0.46, respectively, and are both significant at the 0.05 level. For context, some of the strongest relationships found for teleconnections in this region are presented in Coleman and Rogers (2003), which found (Pearson) correlations between the winter PNA and precipitation in the Ohio River Valley that ranged between −0.55 and −0.7. Our results suggest that Great Lakes winter precipitation is significantly influenced (albeit at a weaker level) by this same teleconnection pattern, which in turn is partially controlled by ENSO and its midlatitude atmospheric response. Interestingly, PLSp1 only projects onto the CP variant of ENSO (CC = −0.32), but is independent of the EP variant (CC = −0.02), consistent with the results in Yu and Zou (2013) over a nearby region that extended further toward the Gulf Coast. The other three PLSp patterns are only weakly and insignificantly correlated to all of the teleconnection indices. Thus, it appears that other major winter precipitation patterns over the Great Lakes basin are related to anomalies in atmospheric flow that are mostly unrelated to major modes of atmospheric variability in the Northern Hemisphere. These unique atmospheric patterns are illustrated by compositing actual and anomalous 500-hPa geopotential heights, 300-hPa wind speeds, and North American precipitation for both the positive and negative phases of each PLSp pattern (Figs. 3–6).
All of the composites show the importance of the position of the jet to Great Lakes precipitation, particularly in the jet entrance and exit regions over North America. Consistent with the correlations in Table 1, the composites for PLSp1 highlight a PNA-like pattern, with a ridge south of Alaska, a trough over western North America, and another weak ridge over the Southeast during positive PLSp1 events (Fig. 3a). These anomalies are associated with a northward displacement of the jet entrance over North America, with anomalously high precipitation across the northern United States and southern Canada and dry conditions in Mexico and the southern United States. During negative PLSp1 events, the wave train is shifted farther to the west and its sign reversed, and the anomaly over the North Pacific is weaker than under the positive phase (Fig. 3c). There are substantial dry anomalies centered in the southern Great Lakes region, but also in the southwest United States. These patterns are similar to composites developed for the strongest El Niño and La Niña years (as defined by the 11 largest and smallest values of DJF Niño-3.4; Figs. 3b,d), but with some key differences, particularly between the negative phase of PLSp1 and strong El Niño events. During the strongest El Niño events, the trough south of Alaska is deeper and farther east compared to negative PLSp1 events and the subtropical Pacific jet extends farther eastward over North America, bringing substantial rainfall to the southern United States and Mexico. The polar jet is also shifted farther north, and negative precipitation anomalies over the Great Lakes are relatively mild. The pressure patterns and precipitation anomalies for PLSp1 are more consistent with a positive PNA, which is associated more with weak-to-moderate versus strong El Niño conditions (Zhang et al. 1996; Hoerling et al. 2001). This suggests nonlinearity in the relationship between Great Lakes winter precipitation and ENSO conditions, with dry conditions during weak and moderate El Niño events that may become less intense under strong El Niño conditions, at least as defined using the Niño-3.4 index. This is consistent with a recent assessment of peak Lake Ontario water levels and ENSO conditions (Carter and Steinschneider 2018). A similar argument for nonlinearity has also been made with respect to Great Lakes winter temperature and ice cover (Rodionov and Assel 2003), but in that case warm anomalies were found to be more consistent under strong El Niño conditions and became more uncertain as event strength weakens.
Composites for PLSp2 show that an east–west dipole in Great Lakes winter precipitation is driven by very different geopotential height patterns and associated shifts in the jet. During negative PLSp2 events, a ridge south of Alaska and trough over western North America cause a northward shift in the jet entrance region and wet anomalies that extend across the Pacific Northwest to Lake Superior (Fig. 4b). However, the elongated trough extends across the entire U.S.–Canadian border and farther east into the Atlantic, with a ridge over southeast Greenland, resembling a negative NAO pattern. The trough center over the Great Lakes drives the core of the jet exit farther to the east, leading to dry conditions in parts of the eastern United States. During positive PLS2 events, the most noticeable feature is a weak ridge over the Hudson Bay, drawing the polar jet farther north (Fig. 4a). This results in storm tracks that are steered up the east coast, causing wet anomalies over Lakes Erie and Ontario and dry anomalies farther to the west over Lake Superior.
The third and fourth PLSp patterns reflect more localized shifts in the jet and regional precipitation (Figs. 5, 6). In both phases of these two PLSp patterns, the jet entrance over North America is shifted to the north, bringing anomalously high precipitation to the Pacific Northwest and generally suppressing precipitation in the Southwest, similar to the negative phase of PLSp2. The major differences between the phases of these two sets of patterns revolve around the position of the jet exit over North America, which determines the location and trajectory of storm tracks over the eastern United States. In general, these tracks are shifted either farther inland (toward Lake Michigan–Huron) or along the coast (toward Lake Ontario), and are influenced by atmospheric anomalies in the Atlantic–North American sector that are asymmetric across different phases of these two PLSp patterns. For the negative phase of PLSp3, the trough over the Hudson Bay and intrusion of cold arctic air blocks storm tracks over the Michigan–Huron basin and drives precipitation farther along the coast, as compared to the positive phase. A similar effect is seen during the negative phase of PLSp4, although the trough is located farther south over the eastern United States, leading to more intense precipitation directly along the coastline and driving the dipole between Lakes Ontario and Erie seen in Fig. 2.
For comparison, we also examined patterns identified through PCA, which are similar to those for PLS, but with some key differences (Fig. S1 in the online supplemental material). For example, the loading vector for the first principal component represents a lakewide average of precipitation, similar to PLSp1, but puts more weight on Lake Ontario and less on Lake Superior. This suggests that precipitation exhibits higher covariance over the middle and lower basin, but the effects of a single mode of atmospheric variability (ENSO) are most consistent across the upper and middle lakes. The second patterns for PLS and PCA are near identical, but the third pattern for PCA consists of alternating positive and negative weights across each of the lakes that are more likely a result of orthogonality constraints than a physically plausible pattern of variability, and are not very consistent with composites of geopotential heights and North American precipitation for the third principal component (Fig. S2). Overall, precipitation patterns defined by the PLS patterns appear more physically consistent with major modes of atmospheric flow that deliver moisture to specific regions of the Great Lakes than the PCA-based counterparts, supporting the use of PLS in this work.
b. Lagged relationship between PLS patterns and sea surface temperatures
To assess the seasonal predictability of the PLSp patterns, Fig. 7 shows composites of SON SSTAs prior to positive and negative phases of each PLSp mode. The composites for PLSp1 show a clear asymmetric response in the tropical Pacific, with strong La Niña conditions in SON preceding positive PLSp1 events, but no obvious El Niño conditions preceding negative events. This is consistent with the results in Fig. 3 and suggests a nonlinear relationship between ENSO conditions and Great Lakes precipitation. SSTA composites for the other PLSp modes show much less consistency, indicating that these modes of variability are more difficult to forecast from SSTs at seasonal lead times. There are, however, a few notable exceptions. In particular, significant SSTAs in the extratropical Atlantic are associated with PLSp2, with warm anomalies near the United Kingdom during positive events and cold anomalies southeast of Newfoundland during negative events. These SSTA patterns reflect those identified in Gastineau and Frankignoul (2015, their Fig. 6), who suggested that these North Atlantic SST patterns in the autumn feed back onto the atmosphere in the early winter to influence circulation over the Atlantic sector. This feedback could form the basis for additional seasonal predictability of Great Lakes winter precipitation, particularly as it relates to the east–west gradient of precipitation across the lakes. In addition, there are moderately significant SSTAs during negative PLSp2 events in the northern Pacific that could also provide some additional seasonal predictability for Great Lakes precipitation. The utility of these signals for out-of-sample prediction will be examined later below.
c. A detailed examination of the relationship between PLS1 and ENSO
The PLSp1 pattern accounts for the largest fraction of variability of Great Lakes winter precipitation (65.1%) and is also the pattern that relates most directly to ENSO. Figure 8 shows rank correlations between DJF PLSp1 and monthly fields of SSTs, OLR, and 500-hPa geopotential height from September to February. While SSTs and heights are available for the entire 1951–2017 period, the OLR data only extend back to the winter of 1975, favoring the use of correlations over composites. We do provide corresponding composite maps in Figs. S3 and S4. Figure 8 depicts the evolution of the planetary-scale wave chain and the underlying oceanic and atmospheric heating that force the atmospheric response, as they relate to PLSp1. In general, an ENSO-driven Rossby wave is clearly discernable, with a PNA-like atmospheric pattern that emerges late in the winter in response to late autumn tropical forcing, with feedbacks to extratropical Pacific SSTs through the “atmospheric bridge” effect (Alexander et al. 2002). However, correlations to tropical SSTs are somewhat weak in autumn because of an asymmetric relationship to SSTs, with strong La Niña conditions during the highest PLSp1 years but little SST response when PLSp1 is lowest (Figs. S3, S4). When examining tropical convection more directly, we find that OLR in a region of the western tropical Pacific is significantly and positively correlated with PLSp1, particularly in October. This relationship is somewhat stronger than that for SSTs, although the time period of analysis is different for the two variables. Unlike the SSTs, there is little evidence to suggest that the strength of the relationship to OLR is asymmetric (Figs. S3, S4). Importantly, although OLR is directly influenced by underlying SSTs, the process of convection is not entirely determined by SSTs, leading to some discrepancies between the two that could be important for prediction. For instance, in the winters of 1980 and 1994, negative OLR anomalies in the western tropical Pacific emerge in the autumn and precede an atmospheric response similar to (although weaker than) that during the strongest negative PLS1 events, but warm SSTs are not prominent in the equatorial Pacific (Fig. S5). Finally, there is also a noteworthy relationship between DJF PLSp1 and October OLR in a region over eastern North America, which coincides with a trough over that same region that is part of a larger wave train over the Pacific–North American region.
The results above show that the ENSO-related wave chain evolves over several months, providing the potential to forecast winter precipitation in the Great Lakes in the autumn based on aspects of ENSO that manifest early in the event life cycle, including equatorial Pacific SSTs and OLR in the tropics. There may also be some predictability from OLR in the midlatitudes over North America, which is surprising given that most seasonal prediction is based on tropical, not midlatitude, forcing. To explore this further, we fit a linear model and a GAM between PLSp1 and October area-averaged OLR in a region between 167.5°E and 177.5°W and between 2.5°S and 2.5°N (OLR.Pa; Fig. 9a), October Niño-3.4 (Fig. 9b), an October ENSO CP index (Fig. 9c), and October area-averaged OLR in a region between 80° and 60°W and between 35° and 60°N (OLR.At; Fig. 9d). We focus on October predictors because this lead time aligns well with the needs of water level management in the Great Lakes system (Carter and Steinschneider 2018). Because the Niño-3.4 and CP indices extend further back than the OLR data, we fit the models for these indices twice based on 1975–2017 and 1951–2017 data. The curve of the fitted GAM for Niño-3.4 exhibits statistically significant nonlinearity (p = 0.05 for 1951–2017 and p = 0.07 for 1975–2017), with larger impacts on PLSp1 under moderate La Niña and El Niño events that become smaller as ENSO conditions become more extreme. By contrast, the GAMs fit using the CP and two OLR indices show no significant nonlinearity. The scatter between PLSp1 and CP is not substantively tighter than that for Niño-3.4, indicating that a focus on central Pacific SSTs independent of the eastern Pacific signal does not resolve the noise with Great Lakes precipitation. Furthermore, accounting for the change in Niño-3.4 between September and November in addition to the October Niño-3.4 value does not improve predictive skill (not shown). Rather, more accuracy is achieved using OLR over the tropical Pacific and eastern North America.
Figure 9e shows leave-one-out predictions of PLSp1 from a series of models over the 1975–2017 time period, including a GAM on Niño-3.4 and linear models for the CP ENSO index, tropical Pacific OLR alone, and OLR from both the tropical Pacific and eastern North America. Under cross validation, the linear model on tropical Pacific OLR outperforms models based on Niño-3.4 and CP ENSO, providing a 9.1% and 8.3% improvement in cross-validated RMSE, respectively. Some of the largest improvements occur in dry winters over the Great Lakes, in particular in 1980, 1994, 1995, and 2003. In 1980 and 1994, OLR was anomalously low in the western tropical Pacific despite neutral SST conditions, as discussed earlier (Fig. S5). In 1995 and 2003, warm SST anomalies were observed, but farther toward the central Pacific and less so in the east, reducing the value of the Niño-3.4 index. The CP index does not predict these years well either, though, because the eastern Pacific is still somewhat warm, while CP ENSO events tend to exhibit a dipole between the west and east tropical Pacific SSTs (Capotondi et al. 2015). The tropical OLR.Pa model also better captures high precipitation years in 1975 and 1976. In 1975, cold SSTs were only modestly below average in October, but OLR anomalies were large, similar to the situation in 1994 but with signs reversed. In 1976, a strong La Niña was emerging in October, but the nonlinearity of the GAM reduces the effect magnitude under these conditions, leading to underprediction by the Niño-3.4 model. The CP index also misses these events.
Interestingly, the model that uses both tropical and North American OLR performs the best under cross validation of all the models tested. It is better able to capture some of the wettest winters, including 1985, 2005, 2006, and 2008, although it still underestimates the most extreme winter events. If the OLR.Pa and OLR.At indices are used in a logistic regression model to predict the sign of PLSp1 anomalies (not shown), out-of-sample model predictions are 73% accurate, suggesting that these predictors could be the basis of a strong categorical forecast. When residuals from the linear regression model based only on OLR.Pa are correlated to the October OLR field, the OLR.At region remains prominent (Fig. S6a). A similar regression against 500-hPa geopotential height shows that OLR.At is associated with a ridge over that same region in October (Fig. S6b), also seen in Fig. 8. However, the OLR.Pa and OLR. At signals are not independent (CC = 0.41). Composites of October 500-hPa geopotential heights and OLR during years with the largest and smallest values for the October OLR.At index show a globally symmetric pattern centered around strong anomalies in the OLR.Pa region (Fig. S7). These OLR and pressure composites show perturbations in the Walker and Hadley circulation cells consistent with a global ENSO response, suggesting that the OLR.At index is capturing an extratropical manifestation of the ENSO process near the Great Lakes. Thus, OLR measures outside of the tropical Pacific may provide the ability to better measure the global reach of ENSO teleconnections prior to the winter season, improving the capacity for seasonal prediction. In fact, similar improvements in out-of-sample prediction were achieved using October OLR measures in the extratropical Southern Hemisphere (not shown), which are unlikely to relate directly to Great Lakes climate, but rather are proxies for the global strength of the emerging ENSO teleconnection.
d. Seasonal forecasts for each of the Great Lakes and relation to OLR event years
Based on Fig. 9, October OLR values from the tropical Pacific and eastern North America are used in linear regressions for winter precipitation over each of the Great Lakes (Fig. 10, Table 2). We also add regional SON SSTs from four regions in the extratropical Pacific and Atlantic (see Fig. 7) to determine if any additional forecast skill can be achieved beyond the aforementioned OLR signals. Several insights emerge from Fig. 10 and Table 2. First, when only the two OLR regions are used, the regression can explain at least 20% and up to 32% of the variability in winter precipitation across the lakes. This prediction skill at a 2-month lead time is impressive and may provide potential to inform lake-level management. Second, the SST regions appear to add little to no additional predictive power over the OLR indices. The SST.NAt index is the only SST-based predictor with a significant coefficient, and only for Lake Superior. Further, under cross validation, the out-of-sample prediction error increases for all lakes when using SST-based predictors. Finally, the OLR.Pa index is only significant in the regression for Lake Superior, while the OLR.At index is significant for all lakes but Superior. Therefore, it seems that the ENSO-related signal from the tropical Pacific is strongest in the western portion of the Great Lakes, at the terminus of an enhanced storm track along the U.S.–Canadian border. The OLR.At index exerts more influence over precipitation in the center and eastern Great Lakes basin, with the strongest effect centered on Lake Erie. The underlying reason for this response is unclear. Although an initial examination of the OLR.At index suggests it is a manifestation of the ENSO process (Fig. S7), a thorough examination of its unique impact on the eastern lakes cannot be provided here due to space limitations, but should be pursued further in another analysis.
Finally, we explore the validity of the continuous modeling framework based on OLR in Fig. 11, which shows climate composites during OLR and non-OLR years for both El Niño and La Niña signals, as defined in Chiodi and Harrison (2015) and Harrison and Chiodi (2017). For El Niño, OLR events exhibit a much broader and eastward OLR extent compared to non-OLR events, as well as a deeper Aleutian low and ridge situated farther east over North America. The precipitation response over North America, including the Great Lakes, appears spatially similar but slightly stronger for OLR events. For La Niña events, there is more asymmetry. Tropical OLR is somewhat broader under OLR versus non-OLR events, but the major difference is the ridge south of Alaska, which is stronger and situated farther west during OLR La Niña events. Accordingly, the strongest precipitation response during OLR events occurs in the Pacific Northwest, but during non-OLR events, this region is relatively dry and the Great Lakes sees the largest increase in winter precipitation. These results are confirmed based on composites of Great Lakes precipitation in each OLR and non-OLR event category: most Great Lakes exhibit significant increases in precipitation during non-OLR La Niña events but are more inconsistent during OLR events (Table 3). For El Niño events, the drying response in the Great Lakes basin is similar between OLR events and non-OLR events and is nearly statistically significant at the 0.1 level. Overall, there may be important differences between Great Lakes precipitation responses under OLR versus non-OLR events, particularly during La Niña, but not necessarily with the strongest response only associated with OLR events, as argued in Harrison and Chiodi (2017). Further work is needed to ensure this asymmetric response is not an artifact of a limited sample, and if confirmed, determine how best to use this asymmetric response in forecasts for the Great Lakes.
This study contributes a diagnostic assessment of major patterns of winter precipitation across the Great Lakes with the greatest covariability with atmospheric flow, with an emphasis on seasonal prediction. Correlation and composite analyses revealed that the first PLS-based mode of precipitation variability is related to ENSO, while the other modes represent unique anomalies in atmospheric flow that govern precipitation gradients over the basin, with limited seasonal predictability. There was also evidence of an asymmetric and nonlinear response of Great Lakes winter precipitation to SST-based measures of ENSO, with a more distinct linkage to La Niña SST conditions and reductions in response magnitude under the strongest SST-based measures of El Niño and La Niña. In addition, the first PLS mode only projects onto the CP variant of ENSO, independent from the EP variant. An examination of the ENSO-related wave train related to PLSp1 highlighted regions of seasonally lagged OLR that were more directly and linearly related to Great Lakes precipitation than SST-based metrics of ENSO. Several linear and GAM models were built to forecast Great Lakes winter precipitation with SST- and OLR-based indices in the tropics and extratropics, with results suggesting that only using OLR linked to the ENSO process provides the most robust forecasting framework. These models generally exhibited distinct precipitation signals that vary in an east-to-west gradient across the Great Lakes basin, although the spatial resolution of these distinctions remains somewhat unclear. Though we observed a linear relationship between tropical Pacific OLR and Great Lakes precipitation, we examined whether this relationship was being driven by a threshold response to a select number of strong OLR events, as reported by Chiodi and Harrison (2015) and Harrison and Chiodi (2017). We found a nonlinear response of Great Lakes winter precipitation to tropical Pacific OLR and non-OLR ENSO events, but with the key distinction that some of the greatest precipitation responses were linked to non-OLR La Niña events. The results of this work indicate that there is potential to improve seasonal forecasts of Great Lakes winter precipitation by optimizing the use of remotely sensed OLR across the globe to better determine the strength of emerging global teleconnections during ENSO events. More broadly, this insight could contribute toward improved subseasonal to seasonal climate prediction in different regions around the globe that experience climate variations linked to the ENSO process. This effort is left for future work.
We would like to thank three anonymous reviewers for their constructive and insightful feedback that helped to improve this work.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-18-0128.s1.