Winter/Spring Runoff Is Earlier, More Protracted, and Increasing in Volume in the Laurentian Great Lakes Basin

Winter/spring runoff has changed in streams worldwide due to climate change, particularly in temperate areas where winter/spring streamflow depends on snowmelt. Such changes potentially affect receiving waters through altered nutrient loading and mixing patterns. The Laurentian Great Lakes are an important freshwater resource and have experienced a myriad of impacts due to climate change. We analyzed 70 years of stream gauge data in the Great Lakes Basin to test for changes in timing, duration, and amount of winter/spring runoff during the period 1950–2019. We found strong evidence for earlier runoff in each of the Great Lakes except Lake Erie, protracted winter/spring runoff throughout the Great Lakes Basin, and a higher runoff depth during the winter‐spring period over time for all watersheds except Lake Superior. Lake Ontario had the greatest change in the date by which 50% of the Jan–May runoff had been discharged (6 days earlier from 1950 to 2019). For winter/spring runoff duration, the most extreme change was observed in Lake Erie (increase of 19 days), and for runoff depth, the greatest change was in the Lake Huron Basin (increase of 3.3 cm). Results were similar for natural and impacted streams. Our results demonstrate dramatic changes in runoff patterns over the last seven decades in the Great Lakes Basin concomitant with previously published changes in precipitation and snowpack. Shifts toward earlier, more protracted, and more voluminous runoff likely change nutrient loading and mixing patterns that influence primary producers, particularly in the nearshore areas of the Great Lakes.


Introduction
Stream hydrology is changing rapidly all over the world with climate change (Blöschl et al., 2017;Musselman et al., 2017;Stewart, 2009).Consequently, receiving water bodies such as lakes may also be impacted on a global scale due to changes in the seasonality of stream flow and resultant changes in inflow volume, timing, and chemical composition.In temperate zones, changes in precipitation patterns and snowmelt patterns have resulted in earlier spring runoff (Barnett et al., 2005;Blöschl et al., 2017;Stewart et al., 2005).Runoff in many areas is also more protracted and may occur in smaller melt events throughout winter rather than as a single snowmelt event during spring (Musselman et al., 2017).Changes in runoff timing may interact with other concurrent physical changes in lakes including altered stratification patterns (Livingstone, 2003), resulting in drastic changes to lake physical conditions during the winter and spring, such as mixing patterns and streamflow intrusion.Other physical changes are also happening concurrently.For example, seasonal ice cover in temperate lakes has decreased in duration (Hodgkins et al., 2002;Sharma et al., 2019) and surface temperatures have increased (O'Reilly et al., 2015;Woolway et al., 2020) while stratified periods have lengthened (Livingstone, 2003).Trends in deepwater temperatures are changing but vary by lake (Pilla et al., 2020).Linking changes in winter/spring runoff timing, duration, and amount to climate change is important to contextualize these other well-known effects of climate change.
Winter/spring runoff, including the snowmelt pulse in temperate lakes, can be an influential driver of lake processes throughout the year.The snowmelt pulse often brings high nutrient loads that drive growth of primary producers (Isles et al., 2017;Michalak et al., 2013;Rosenberg & Schroth, 2017).In some systems, winter/spring runoff may also be a significant part of the yearly water budget (Mielko & Woo, 2006) and the amount of runoff may vary widely from year-to-year depending on snowpack (Sadro et al., 2018).Furthermore, the timing and temperature of water inputs in relation to ice cover and thermal stratification may dictate how inputs are mixed into the water column when they enter the lake.For example, runoff that enters a lake during ice cover and cryostratified conditions (i.e., winter stratification - Yang et al., 2021) would likely remain as a plume directly under the ice where it is more likely to travel toward an outflow (Cortés et al., 2017).In contrast, snowmelt runoff entering during ice free and isothermal conditions or ice free and stratified conditions (Roberts et al., 2018) is much more likely to contribute nutrients that are retained in the lake and affect lake biogeochemistry.Therefore, the timing of the bulk of winter/spring runoff has strong implications for the fate of associated nutrients once they enter the lake.
The Laurentian Great Lakes (hereafter Great Lakes) are in a region with extensive snow accumulation and substantial winter/spring runoff, but their physical dynamics have been changing in recent decades.Ice cover has decreased in the Great Lakes Basin since at least the 1980s (Baijnath-Rodino et al., 2018), while snowfall trends have been more difficult to interpret.Snowfall decreased from 1980 to 2015 north of lakes Superior and Huron but increased in many areas further south (Baijnath-Rodino et al., 2018).However, earlier historical data indicated that snowfall increased around lakes Superior, Michigan, and Huron, but not Erie and Ontario (Kunkel et al., 2009).The average snow depth that accumulates and remains during winter in the Great Lakes Basin has decreased over time in concert with increases in air temperature (Suriano et al., 2019).The decrease in snowpack, despite increases in snowfall in some areas, implies that some snow accumulation is melting and entering lakes as runoff earlier in the year rather than accumulating as snowpack.
Decreased ice cover basin-wide, changes in winter precipitation toward more rain and less snow, and potentially increased lake-effect snow in the most northern parts of the basin due to decreased ice cover are expected in the future (Notaro et al., 2015).Simulations predict increased winter precipitation overall with a higher proportion of winter precipitation as rain on the leeward sides of the Great Lakes (Notaro et al., 2015).Changing winter/spring runoff timing may have major impacts on productivity in the Great Lakes.One example of this is Western Lake Erie, where the nutrient dynamics of winter/spring runoff have been studied particularly well; spring floods deliver high nutrient loads that contribute to harmful algal blooms during the summer (Michalak et al., 2013;Stumpf et al., 2012).Therefore, further elucidation of how climate change and its associated changes in snowpack and precipitation have affected winter/spring runoff patterns are likely to identify processes that need further study to improve the interoperation and simulations of future changes in the Great Lakes.
We examined how runoff has changed over time in the Great Lakes Basin.Analyses of historical changes in runoff have been limited in their spatial range in the Great Lakes Basin (e.g., Byun et al., 2019;USGS data only).Other runoff studies have mainly focused on future projections of hydrological conditions in the Great Lakes basin or surrounding areas (Byun et al., 2019;Cherkauer & Sinha, 2010;Demaria et al., 2016;Naz et al., 2016).To our knowledge, no previous studies have quantified changes in winter/spring runoff timing, duration, and amount for the entire Great Lakes Basin over several decades.We tested three hypotheses related to winter/spring runoff for each of the Great Lakes from 1950 to 2019: (a) the timing of winter/spring runoff is earlier than past decades, (b) the duration of runoff is more protracted than past decades, but (c) the amount of runoff has not changed over time.These hypotheses were tested using historical stream gauge data from the US and Canada.

Materials and Methods
The Great Lakes are the largest chain of freshwater lakes in the world and are in a region with extensive snow accumulation and substantial winter/spring runoff.Each of the Great Lakes may vary in their response to climactic conditions due to unique combinations of bathymetry, latitude, land use, and hydrological conditions.Residence times vary from approximately 2.7 years in Lake Erie to 173 years in Lake Superior (Quinn, 1992).
Agricultural activity is greater in the southern parts of the Great Lakes Basin than in the northern parts (Lunetta et al., 2010), and the Lake Erie and Lake Michigan sub-basins are the most populated (Méthot et al., 2015).
Stream gauge data for the entire Great Lakes Basin were collected from US Geological Survey (USGS) and Environment and Climate Change Canada (ECCC) public data.We downloaded stream gauge data for all stream gauges, organized by Hydrologic Unit Codes (HUC 8 level in USA) or watersheds (tertiary watershed level in Canada) in the Great Lakes Basin following a published list of HUCs and watersheds (Table 4 in Neff et al., 2005; see also Meyers et al., 2023) from the inception of gauging programs through 2019.The USGS HUC 8s and ECCC tertiary watersheds are hereafter referred to as "sub-watersheds" and they include several stream gauges each.Data were extracted from the ECCC Historical Hydrometric Data web site (https:// wateroffice.ec.gc.ca/mainmenu/historical_data_index_e.html) on 12 September 2020 and from the USGS Daily Values Web Service (https:// waterservices.usgs.gov/rest/DV-Test-Tool.html) on 27 September 2020.
Data were cleaned to retain only years and stream gauges that had sufficient data for analyses.We retained only data from 1950 to 2019 because more stream gauges were operational after 1950 (Figure 1) and because Neff et al. (2005) used a similar cutoff date of 1948 to calculate historical base flow in the Great Lakes Basin.Rivers that drained to Lake St. Clair were included in the Lake Erie watershed.We eliminated sites that did not have watershed area included in their metadata.Rivers that connect two great lakes (St.Mary's, Detroit, St. Clair, and Niagara Rivers) were excluded from the analysis because we were more interested in watershed inputs than the flow of water from one Great Lake to another.We eliminated provisional data from USGS, which was predominantly unverified recent data at the time of download (late 2019).We retained data when the gauged rivers were ice covered because, although agencies handle ice cover differently, they are handled in the same manner at any single site.For example, stream gauges in some US states have data removed during ice cover, while data from gauges in other states are estimated during ice cover.Days with no data were excluded from analyses.Data from all gauges were included if there were multiple stream gauges in one sub-watershed because our indices primarily dealt with runoff timing rather than volume, and the index derived from volume (runoff depth) took gauged area into account.Although a small number of streams have multiple gauges, the effects on results should be negligible and have been combined this way in previous research (Neff et al., 2005).We used data from a total of 1501 stream gauges.
Runoff indices (timing, duration, and amount) were calculated for all stream gauges for all years.We defined the winter/spring window as January 1-May 31.Although seasonality may be different among sub-basins and years, this window largely captures the winter spring runoff entering temperate lakes including any spring snowmelt and is therefore a robust description of the winter/spring period (Hrycik et al., 2021).To calculate runoff timing, we used the winter/spring center of volume (WSCV), which measures the day of year (as Julian Day) when half of the total winter/spring discharge has occurred (Burns et al., 2007;Zion et al., 2011).To quantify how protracted runoff was in a given year and for a given stream gauge, that is, the duration of winter/spring runoff, we calculated the inter-quartile distance (IQD).The IQD is defined as the number of days that pass between the dates at which 25% and 75% of the winter/spring discharge has occurred (Hrycik et al., 2021).Runoff depth during the winter/ spring period was calculated as the cumulative volume of water (m 3 ) measured at a stream gauge during the winter/spring period divided by total area of the watershed associated with each stream gauge, which enables comparisons between gauges from watersheds with different areas.This also enabled comparisons within tributaries with multiple gaged sub-basins.Runoff depth, therefore, represents the depth of water lost from the entire watershed area from January through May of a given year.Each of these indices (runoff timing, duration, and depth) was calculated for every stream gauge in every year.Runoff depth data had some extremely high outliers in specific gauges which altered trends.These outliers were orders of magnitude higher than other data and likely represented unusual circumstances at the stream gauges rather than actual high water years.Consequently, we removed years and sub-watersheds with runoff depth greater than the 99.5th percentile of all data from analyses.
After the three runoff indices were calculated for each stream gauge and year, we calculated weighted average indices for each sub-watershed (HUC 8 or tertiary watershed) by year.We did not examine changes over time for each stream gauge because most gauges did not record for the entire period of study, so a sub-watershed average for each year gave more complete results.When calculating averages runoff timing and duration indices for each stream gauge were weighted by their watershed area so that gauges that captured runoff from a greater area had more influence than gauges on minor streams.Runoff depth was directly averaged because by definition it was already normalized by area.To justify combining stream gauges in this manner, we examined means and ranges of runoff indices from single stream gauges compared to sub-watershed indices to make sure results were similar among stream gauges.
Responses in runoff indices over time were addressed using linear regression.We fitted linear regressions for each runoff index (WSCV, IQD, and runoff depth) over years for each sub-watershed.We then calculated mean slopes from all the sub-watersheds within each lake basin to obtain one average value for each Great Lake and calculated bootstrapped 95% confidence intervals around the mean slope by resampling the calculated slopes with replacement 1000 times (Hrycik et al., 2021).The slope estimates from each sub-watershed were then used to calculate the average change in each runoff index over time for each of the Great Lakes.Statistical significance for change in indices over time was determined by whether the bootstrapped confidence intervals overlapped a slope of zero, that is, confidence intervals that did not overlap zero were considered significant.
We analyzed January-May precipitation data to understand whether precipitation trends were likely to drive trends in runoff depth over time.We used average monthly overland precipitation estimates described in Hunter et al. (2015) (downloaded from https://www.glerl.noaa.gov/pubs/tech_reports/glerl-083/UpdatedFiles/ on 21 November 2023).We summed monthly estimates for each Great Lake watershed (reported in mm) from January through May to derive total precipitation estimates, and then ran linear regressions over time on January-May precipitation for the years in our study .
We evaluated the effects of human impacts on runoff indices by separating stream gauge data into "natural" and "impacted" categories based on available USGS and ECCC metadata."Natural" streams had less regulation and lower human impact, and were designated as reference streams by USGS (Mai et al., 2022) and natural streams ECCC in gauge metadata."Impacted" streams were those designated as non-reference streams by USGS or regulated streams by ECCC.Within natural and impacted categories, we re-calculated all runoff indices and repeated bootstrapping analyses for testing significance.Streams that were not categorized in USGS or ECCC metadata (10.8% of data) were excluded from this analysis.

Results
A comparison of indices from single stream gauges within each sub-watershed justified averaging runoff indices from multiple stream gauges within a sub-watershed.Runoff timing differed by stream within each subwatershed, but WSCV and IQD calculated from different stream gauges within the same sub-watersheds were similar, supporting our choice to calculate averages of stream gauges in a sub-watershed to get yearly runoff metrics.For example, the annual average range for the WSCV for all the gauges within each sub-watershed was 10 days, while the IQD was slightly more variable (average range within year and sub-watershed = 15 days).However, the medians and means of both WSCV and IQD were extremely close (average difference between median and mean for WSCV = 1.0 days and for IQD = 1.4 days), indicating that distributions were not heavily skewed or influenced by outliers.Therefore, we are confident to represent runoff indices for each year and subwatershed as an average of the indices from all stream gauges (weighted by gauged area) within the subwatershed.
The timing of winter/spring runoff was dependent on the Great Lake watershed and its sub-watersheds, but tended to become earlier over time.The bulk of runoff happened earliest in the year in lakes Ontario and Erie as indicated by the WSCV, while Lake Superior had the latest runoff in the year and lakes Huron and Michigan were intermediate (Figure 2).Most sub-watersheds of all Great Lake basins experienced a trend toward earlier runoff, but this change toward earlier runoff was significant only in lakes Ontario and Michigan.In Lake Ontario, the WSCV advanced to 6.0 days earlier and in Lake Michigan WSCV advanced 4.3 days earlier over the seven decades of the study.However, even though only lakes Michigan and Ontario showed significant trends in the basin-wide mean slope, the majority of sub-watersheds had negative slopes for the relationship of WSCV over time, except in the Lake Erie drainage (Figure 3; Table 1).Notably, 82% of the sub-watersheds in the Lake Superior drainage and 74% of the sub-watersheds in the Lake Huron drainage showed a trend toward earlier runoff, even though the mean basin slopes were not significantly different from zero (Table 1).Runoff timing changes were weakest in Lake Erie, which had a non-significant 95% confidence interval and only 59% of sub-watersheds trended toward earlier runoff.WSCV over time varied spatially, with sub-watersheds near western Lake Erie showing a particularly high density of positive slopes, indicating that runoff is getting later near western Lake Erie in contrast to many other areas (Figure 4).
Winter/spring runoff entering the Great Lakes trended to be more protracted over time.IQD was similar among the Great Lakes watersheds and most sub-watersheds trended toward higher IQD, with notable exceptions in some sub-watersheds in Superior, Huron, and Ontario (Figure 5).The mean basin trend in IQD was significantly positive in Lakes Michigan, Ontario, and Erie, with 97%, 90%, and 100% of respective sub-watersheds showing a trend toward more protracted runoff (Table 1; Figure 3).IQD increased by 9.4 days in Lake Michigan, 19.4 days in Lake Erie, and 14.9 days in Lake Ontario from 1950 to 2019.Although the mean slope was not significantly different from zero in Lakes Superior and Huron, sub-watersheds for each lake were predominantly trending positive (74% of sub-watersheds for Lake Huron and 77% for Lake Superior; Table 1).Most exceptions that showed decreased IQD were found near northern Lake Huron (Figure 4).
On average, the winter/spring runoff depth increased over time (Table 1, Figures 3 and 6).Basin average increases in runoff depth were significant in lakes Michigan, Erie, and Huron.Winter/spring runoff depth increased by 3.3 cm in the Lake Huron Basin, 2.5 cm in the Lake Michigan Basin, and 3.5 cm in the Lake Erie Basin from 1950 to 2019.Although the average slope was not significant, most sub-watersheds in Lake Ontario also experienced increased runoff volume (71% of sub-watersheds).Lake Superior did not exhibit a significant trend, with 59% of sub-watersheds experiencing trends toward decreasing runoff, while 41% of sub-watersheds trended toward increased runoff (Table 1; Figure 4).January-May precipitation increased over time in most of the Great Lakes Basin, including significant increases from 1950 to 2019 for Lakes Huron (linear regression, F 1,68 = 7.67; p = 0.01), Michigan (F 1,68 = 8.11; p = 0.01), Erie (F 1,68 = 4.13; p = 0.05), and Ontario (F 1,68 = 3.96; p = 0.01).Lake Superior was the only watershed that did not experience an increase in January-May precipitation over the period of our study (F 1,68 = 0.46; p = 0.50).
Impacted and natural streams had similar trends in runoff indices, although significance levels varied (Table 1).Some trends became non-significant for one or both datasets when the larger dataset was split into impacted and natural sub-watersheds (i.e., WSCV for lakes Michigan and Ontario, IQD for lakes Michigan and Erie, and runoff depth for lakes Huron, Michigan, and Erie).Others had significant sub-sets of data when the full data set was nonsignificant (i.e., WSCV in Lake Erie, IQD in Lake Huron, and runoff depth in Lake Ontario).Nonetheless, almost every grouping had the same trends when split into impacted and natural sub-watersheds, as evidenced by the percentages of positive and negative slopes for trends in runoff indices over time (Table 1).The only runoff index that consistently showed changed trends between impacted and natural sub-basins was runoff depth for Lake Michigan.Seventy-one percent of all sub-watersheds for Lake Michigan had a positive trend toward increasing runoff over time, while only 45% of sub-watersheds with natural streams had increasing trends.Overall, trends  were strikingly similar in all groups, with every grouping showing a majority of negative trends for WSCV and positive trends for IQD, indicating earlier and more protracted runoff over time that was largely independent of basin development.Runoff depth was slightly more variable, with all but two of the possible groupings showing a positive trend toward increased runoff over time (Table 1).The sample size for impacted streams was larger than that for natural streams for all of the Great Lakes except Lake Superior.The number of impacted and natural subwatersheds listed is not equivalent to the total number of sub-watersheds because some of the watersheds contained gauges on both impacted and natural streams.

Discussion and Conclusions
Overall, winter/spring runoff became earlier, more protracted, and had higher volume in the Great Lakes Basin over the period 1950-2019.Runoff timing was significantly earlier in Lakes Michigan and Ontario, and nonsignificant in Lakes Superior and Huron, despite strong trends toward becoming earlier.Lake Erie was the only Great Lake that showed only weak evidence of earlier runoff.Our results indicate that the fastest change toward earlier runoff occurred in the Lake Ontario watershed.The duration of winter-spring runoff became significantly longer in Lakes Erie, Ontario, and Michigan, and similar but non-significant trends were found in Lakes Superior and Huron.The amount of winter-spring runoff increased in all the lakes except Lake Superior.Increased runoff was significant in Lakes Huron, Michigan, and Erie, but also showed an increasing trend in Lake Ontario.Changes in runoff have implications for mixing patterns, nutrient dynamics, water balance, and potentially nearshore primary production in the Great Lakes.
Little difference was evident between impacted and natural streams, indicating that our results are likely due to broad changes associated with climate change rather than local alterations to sub-watersheds.Several land use changes have the potential to alter runoff or stream gauge readings, including irrigation (Chang et al., 2015;Destouni & Prieto, 2018), drain tiling (Sloan et al., 2016), and changes in vegetation (e.g., removing vegetative material; Brooks et al., 1994).Changes within the streams themselves also may impact gauge readings and water movement, including channelization (Anim et al., 2018) and regulation of dams (Zimmerman et al., 2010).Land use changes and in-stream changes may drive some outliers in our dataset, however, they do not seem to be driving overall trends based because responses in natural and impacted streams were similar.Climate change Note.Percent negative/positive refers to the percentage of watersheds that showed a negative or positive slope, mean slope is the average slope across all sub-watersheds for the response variable over years for the period 1950-2019, and 95% CI is the bootstrapped 95% confidence interval around the mean slope.Values for all sub-watersheds in a lake basin (impacted plus natural streams) are bolded.Significant slopes (95% CIs that do not overlap zero) are marked with *.The number of impacted and natural watersheds does not add up to the total number of watersheds because some watersheds have both impacted and natural streams, and some streams were not classified.

Water Resources Research
10.1029/2023WR035773 HRYCIK ET AL. 8 of 14 operates on a global scale and has been implicated in changes in runoff timing and precipitation patterns in many regions (Musselman et al., 2017(Musselman et al., , 2018) ) and this is likely the case in the Great Lakes basin as well.
Runoff dynamics, including runoff timing, duration, and amount, drive allochthonous nutrient loading in lake systems (Isles et al., 2017;Michalak et al., 2013;Rosenberg & Schroth, 2017).External nutrient loading from runoff has been implicated in worsening algal blooms in Lake Erie (Michalak et al., 2013;Stumpf et al., 2012), but is also important for phytoplankton growth in the other Great Lakes (e.g., Johengen et al., 2008) and many smaller systems (Lathrop, 2007;Xu et al., 2010).Winter/spring runoff often brings a high load of bioavailable nutrients (Rosenberg & Schroth, 2017;Seybold et al., 2022;Sickman et al., 2003) that can fuel primary production (Michalak et al., 2013;Slemmons & Saros, 2012).Effects of changes in runoff dynamics in the Great Lakes are likely most noticeable in nearshore areas, where the spring development of a thermal bar limits mixing between nearshore and offshore areas (Rao & Schwab, 2007;Spain et al., 1976) and prevents nutrients in river inflow from mixing into deeper pelagic zones in the Great Lakes (Makarewicz et al., 2012).Phytoplankton biomass tends to be higher nearshore due to the thermal bar (Pavlac et al., 2012;Pothoven & Vanderploeg, 2020) and we expect that more protracted runoff that begins earlier in winter before the thermal bar develops will result in a less dramatic contrast between nearshore and offshore phytoplankton blooms in spring.
The timing and magnitude of snowmelt may also have implications for how inflow impacts flushing rates and mixing in the Great Lake nearshore zones.For example, the Trent River is the main driver of the flushing rate in the Bay of Quinte on Lake Ontario (Shore, 2020).With more protracted runoff and an increased runoff volume, the flushing of the Bay of Quinte and other bays may occur over a longer time period, and also be affected by the presence of ice cover.The interaction of the timing of runoff entering a lake and the lake's antecedent conditions are critical for the fate of the external nutrient load, particularly in nearshore areas that are ice-covered in winter.Cold runoff that occurs during ice covered, cryostratified periods (Yang et al., 2021) during winter may travel as a plume at the surface of the lake just under the ice (Bergmann & Welch, 1985;Cortés et al., 2017).When ice cover is not present, early snowmelt that enters lakes during spring mixing may enter as a thick, unconfined surface plume, while late snowmelt following the onset of thermal stratification may come in as a thinner, plunging plume (Roberts et al., 2018).Therefore, earlier runoff events during winter ice cover may travel under the ice as a plume during winter inverse stratification or be mixed into the water column when the runoff occurs during isothermal spring mixing, as we often expect with a "typical" spring flood.We expect that earlier winter runoff pulses such as small snowmelt events that enter ice-covered bays would travel further out into the lake as a plume if winter stratification is present (Cortés et al., 2017) or become mixed into the water column if under-ice convection is taking place just prior to ice-out (Bruesewitz et al., 2015).
These differences in the flow paths of runoff, depending on hydrothermal condition of the lake, may impact phytoplankton growth later in the year if the uptake or sequestration of nutrient is altered.Smaller snowmelt pulses that occur under ice are expected to decrease summer phytoplankton biomass in smaller temperate lakes (Hrycik et al., 2021;Pierson et al., 2013) and the same may be true for nearshore areas or small bays in the Great Lakes.However, future increases in the duration and intensity of thermal stratification and warmer temperatures during summer are expected to intensify phytoplankton production (O'Neil et al., 2012;Paerl & Huisman, 2008).Contrasting expectations of the role of early-season runoff entering ice covered lakes versus that occurring during the open-water season make interpretation and prediction of primary production trajectories with climate change difficult, and the effects of runoff timing on lake productivity may also differ by region.Many of the changes we observed in runoff timing and volume result from changes in precipitation and snow storage patterns in the Great Lakes Basin (Suriano et al., 2019;Zhang et al., 2019).Earlier and more protracted runoff can be explained by episodes of snowmelt throughout the winter, particularly as a consequence of rain-onsnow events, while increased runoff volume is likely the result of increased precipitation as rain and runoff occurring at times when evapotranspiration is low early in the season (Kiefer et al., 2019).Our precipitation analysis indicated that January-May precipitation increased through the period of our study, and other research indicates that this trend will continue into the future with climate change.Precipitation in winter and spring is expected to increase with climate change in the Great Lakes Basin due to higher rainfall during more extreme events (d'Orgeville et al., 2014;Hayhoe et al., 2010), increased lake-effect precipitation downwind of each of the Great Lakes (Notaro et al., 2015), and increased total winter precipitation compared to current conditions (Zhang et al., 2020).A shift toward earlier runoff is widespread in temperate areas beyond the Great Lakes (Rauscher et al., 2008;Stewart et al., 2005).Magnitude of runoff, however, varies by region (e.g., Stahl et al., 2010).
Runoff timing for Lake Erie became more protracted and Lake Erie had the most dramatic increase in volume of runoff of all the Great Lakes, but the timing of the center of winter/spring runoff did not change as much in Lake Erie as the other lakes.We suspect that the lack of change in runoff timing was due, in part, to a lower importance of snowmelt storage to the hydrological cycle of the Lake Erie watershed, which is at the lowest latitude of all the Great Lakes and has the most limited snowpack (Suriano et al., 2019).Some outliers were evident in our runoff indices, but we chose to broadly keep data in the analysis except for some runoff depth values that were suspiciously high.As a result, runoff indices were extremely variable.Some outliers, but not all, came from shorter data sets where trends were strongly affected by one or a few extreme points.We chose to keep shorter time series despite this variability to have the most complete dataset possible to elucidate broad trends.Furthermore, outlier definitions are extremely variable and sometimes unclear among authors (Aguinis et al., 2013;Sullivan et al., 2021).We did not see step changes or other anomalies that suggested data were in error.Additionally, our analysis of impacted and natural sites did not reveal many differences between streams with high and low human impact, suggesting that trends were consistent among sub-watersheds with different land use changes.
We found that runoff changed dramatically in the Great Lakes Basin over the past 70 years: the majority of winterspring runoff occurs earlier in the year, runoff is more protracted throughout the winter and spring, and runoff is occurring at higher volume over time.Earlier, more protracted runoff due to climate change is common in many regions of the world where snow plays a role in the annual hydrologic cycle (Musselman et al., 2017;Stewart et al., 2004Stewart et al., , 2005)), but higher runoff volumes calculated here reflects local changes in precipitation patterns specific to the Great Lakes Basin (d'Orgeville et al., 2014;Hayhoe et al., 2010).The changes in runoff documented here imply that nutrients are entering the Great Lakes earlier in the season when phytoplankton are potentially less able to use them (Hrycik et al., 2021) and before onset of the spring thermal bar, suggesting declining primary productivity with less well-defined seasonality.Recent oligotrophiction in the Great Lakes has been attributed to reduced nutrient loading and filtering of phytoplankton by introduced Dreissena mussels ( Barbiero et al., 2012;Evans et al., 2011).Our findings suggest that changes in timing, duration, and amount of nutrient loading occurring during the Winter-Spring period may also need to be considered to understand the ongoing changes to Great Lake primary productivity.

Figure 1 .
Figure 1.The number of gauged streams in each Great Lake watershed over time.Black lines indicate US data and gray lines indicate Canadian data.The vertical dashed line indicates 1950, when we began using data for this study.Note that some pre-1900 data are available, but they are sporadic and not shown.

Figure 2 .
Figure 2. Trends in winter/spring center of volume (WSCV), a measure of runoff timing, in the Laurentian Great Lakes from 1950 to 2019.Each color represents a different sub-watershed.Lines are results of linear regressions for each sub-watershed.

Figure 3 .
Figure 3. Distributions of slopes (in days/year) for winter/spring center of volume (WSCV) and inter-quartile distance and runoff depth (m/year) over the period of study within each sub-watershed represented with kernel density estimation, that is, these distributions represent the slopes of the linear regressions represented in Figures 2, 4 and 5. Dashed lines indicate slope of zero (no change) and shaded regions show the bootstrapped 95% confidence interval.Positive slopes indicate later runoff, more protracted runoff, or higher runoff volume over time, while and negative slopes indicate earlier runoff, more condensed runoff, or decreased runoff depth over time.Note that axis scales differ by panel.

Figure 4 .Figure 5 .
Figure 4. Slopes of winter-spring center of volume (WSCV), inter-quartile distance (IQD), and runoff depth over time in units of days/year for WSCV and IQD and m/year for runoff depth.Each point represents a subwatershed.Red points indicate negative slopes, blue points indicate positive slopes, and point size indicates the relative magnitude of the slope.

Figure 6 .
Figure 6.Trends in the runoff depth from 1950 to 2019.Each color represents a different sub-watershed.Lines are results of linear regressions for each sub-watershed.

Table 1
Results of Linear Regressions for Each Runoff Variable, Averaged by Sub-Watershed and Impacted or Natural Streams