HYDROLOGIC VARIATIONS OWING TO SNOWMELT CHANGES IN THE MID LATITUDES By Chanse M. Ford A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degrees of Geological Sciences—Doctor of Philosophy Environmental Science and Policy—Dual Major 2022 ABSTRACT HYDROLOGIC VARIATIONS OWING TO SNOWMELT CHANGES IN THE MID LATITUDES By Chanse M. Ford Snowmelt is a critical hydrologic process in high latitude, non-alpine settings. The water stored in seasonal snowpacks melts in the spring months leading to increased spring streamflow and creating shallow groundwater recharge that helps sustain streams throughout the year from the contributions of baseflow. Many regions across the globe have experienced changes to snowpack dynamics and melt patterns due to increased winter temperatures resulting from global climate change. Currently, most of the research into the changing snowmelt hydrology has been focused on mountainous regions where snowpacks make up larger portions of those regions’ annual water budget. There is little research in these mid-latitude, non-alpine areas and the available research focuses on small areas or examines only one component of the hydrologic system. These understudied regions that receive seasonal snowfall require more thorough examination as changes to winter and spring snow can have negative societal consequences, especially in one of the world’s largest freshwater reservoirs of the Great Lakes. This dissertation contributes to the scientific knowledgebase regarding snowmelt dynamics in non-alpine settings. Novel statistical analyses are utilized to assess the amount of change to winter temperatures and the effects on snowmelt hydrology across spatial scales and decades of observational data. The results from these analyses are then used as a lens to simulate the landscape hydrology to quantify changes to shallow groundwater recharge, which is difficult to assess from empirical data alone. These findings also lead to an examination of the potential economic effects resulting from changes to the snowpack. Chapter 1 lays the groundwork for this dissertation by describing the relevance of the research and gives a brief overview of the different components. The foundational methodology is developed in Chapter 2, where a combination of observed physical data and outputs from several snow and precipitation models are used to classify winters in Michigan from 2003-2017 as warm or cool and quantify the hydrologic changes in those different winter types. The results show warmer winters had less overall snow, which melted earlier contributing to earlier and lower spring stream flows and increased net recharge of groundwater. Chapter 3 then takes this methodology and applies it to the entire eastern portion of the United States that receives seasonal snow from 1960-2019, with results similar to the preceding chapter, and demonstrating that these snow hydrology changes are not limited to Michigan or to the more recent decades. These findings then culminate in Chapter 4, where the Landscape Hydrology Model simulates the snowpack, surface flows and groundwater recharge in Michigan from 2000-2019. These fully distributed simulations show the decreased snow and periodic melting in warm winters has led to increased groundwater recharge and decreased surface flows. Chapter 5 concludes this dissertation by examining the downhill ski industry in Michigan using industry statistics and operational data from the Shanty Creek Resorts, describing the potential economic challenges in a warming future. This work is dedicated to my grandfather, Terry Bishop, who passed away in September 2020 and was only able to read the first chapter. His influence instilled in me the dedication, confidence and curiosity necessary to complete this work and succeed as a scientist. He is sorely missed. iv ACKNOWLEDGMENTS First and foremost, I would like to thank my advisor, David Hyndman. When I originally applied to the doctoral program in the Earth and Environmental Science department at Michigan State University, it was originally to a different lab. Dr. Hyndman contacted me and asked me to consider instead working in the Hydro Lab. Ever since he has been invaluable in his expertise and his support, allowing me to develop this project and make it my own. Without his guidance this body of work would not have been possible. I also would like to thank Anthony Kendall for his vast amount of assistance through the years. Dr. Kendall was heavily involved with every facet of this research, from research ideas to publication. Through his teachings I’ve improved every technical skill necessary to be an environmental scientist. I am a better researcher and science communicator because of his influence. Numerous peers in the department have helped me every step of the way by providing their code, reviews of my work and moral support. Brent Heerspink and Autumn Parish provided support in too many areas to list. Quercus Hamlin, Jake Rousch and Jill Deines were invaluable in honing my coding and geospatial analysis skillset. Many others including Alex Kuhl, Luwen Wan, Jacob Stid, Jeremy Rapp and Sherry Martin were always available as another set of eyes on manuscripts and conference presentations. The love, support and encouragement of my family has always been foundational in my academic success, and the completion of my dissertation is the culmination of that support. My soon to be wife, Martha Timmons, has been here since the beginning of my graduate career as my rock. My parents, Colleen and Dave Wolford, have been supportive of everything I do and v have taught me so much about how to succeed in life. The rest of my family has always been there to help and encourage me in whatever ways I need, including my brother, Bailey Ford, my grandparents, Terry and Kathy Bishop, my aunt, Shannon Braesseale, and my uncle, Sean Bishop. I appreciate all of you more than I can express. This research would not have been possible without the funding I’ve received from a variety of sources of the years. This includes grants from the USDA-NIFA project, the NOAA Tipping Points project and the NASA Invasives Project. Additional funding came from the NOAA-CIGLR Summer Fellows program and Dr. Yao Hu at the University of Delaware for the Runoff Risk Project. Other support in the form of fellowships and teaching assistantships came from the Michigan State University Graduate College, the Department of Earth and Environmental Science and the Environmental Science and Public Policy program. Finally, I would like to acknowledge my undergraduate advisor at the University of Southern Indiana, Dr. Paul Doss. Without his inspiration and guidance, I would not have pursued a career in hydrologic science. vi TABLE OF CONTENTS LIST OF TABLES ……………………………………………………………………………….ix LIST OF FIGURES……………………………………………………………………………….x KEY TO ABBREVIATIONS……………………………………………………………………xii CHAPTER 1: INTRODUCTION ................................................................................................... 1 CHAPTER 2: EFFECTS OF SHIFTING SNOWMELT REGIMES ON THE HYDROLOGY OF NON-ALPINE TEMPERATE LANDSCAPES ............................................................................. 6 Abstract ....................................................................................................................................... 6 1. Introduction ............................................................................................................................. 6 2. Methods ................................................................................................................................ 10 2.1 Study Region................................................................................................................... 10 2.2 Data Sources ................................................................................................................... 13 2.3 Analysis Methods............................................................................................................ 15 2.3.1 Spatial Aggregation ............................................................................................... 15 2.3.2 Multimetric “Warm” and “Cool” Year Classification ........................................... 16 2.3.3 Regional-Scale Snowpack and Streamflow Analyses ........................................... 17 2.3.4 Winter and Spring Net Recharge Analysis ............................................................ 19 2.3.5 Dataset Distribution Tests ...................................................................................... 19 3. Results ................................................................................................................................... 20 3.1 Warm and Cool Year Classification ............................................................................... 20 3.2 Precipitation Analysis ..................................................................................................... 21 3.3 Snowpack Analysis ......................................................................................................... 22 3.4 Streamflow Analysis ....................................................................................................... 26 3.5 Winter and Spring Recharge Analysis ............................................................................ 28 4. Discussion ............................................................................................................................. 31 5. Conclusions ........................................................................................................................... 36 Acknowledgments .................................................................................................................... 37 APPENDIX ................................................................................................................................... 39 CHAPTER 3: SNOWPACKS DECREASE AND STREAMFLOWS SHIFT ACROSS THE EASTERN US AS WINTERS WARM........................................................................................ 42 Abstract ..................................................................................................................................... 42 1. Introduction ........................................................................................................................... 43 2. Methods ................................................................................................................................ 47 2.1 Study Region................................................................................................................... 47 2.2 Data Sources ................................................................................................................... 49 2.3 Analysis Methods............................................................................................................ 50 2.3.1 Spatial and Temporal Aggregation for Analysis ................................................... 50 2.3.2 Multimetric “Warm” and “Cool” Snow Season Classification ............................. 51 vii 2.3.3 Regional-Scale Snowpack, Streamflow and Precipitation Analyses ..................... 52 2.3.4 Distribution Analysis and Significance Tests ........................................................ 53 3. Results ................................................................................................................................... 54 3.1 Winter Type Classification and Trends .......................................................................... 54 3.2 Snow Depth and Precipitation Analysis ......................................................................... 57 3.3 Streamflow Analysis ....................................................................................................... 61 4. Discussion ............................................................................................................................. 65 5. Conclusions ........................................................................................................................... 70 Acknowledgments .................................................................................................................... 71 APPENDIX ................................................................................................................................... 72 CHAPTER 4: SNOWMELT DYNAMICS FROM WARMER WINTERS INCREASES GROUNDWATER RECHARGE IN MICHIGAN ...................................................................... 82 Abstract ..................................................................................................................................... 82 1. Introduction ........................................................................................................................... 82 2. Materials and Methods.......................................................................................................... 85 2.1 Study Region................................................................................................................... 85 2.2 Landscape Hydrology Model (LHM) ............................................................................. 87 2.3 Data Sources ................................................................................................................... 90 2.4 Model Calibration ........................................................................................................... 92 2.5 Analysis Methods............................................................................................................ 93 3. Results ................................................................................................................................... 94 3.1 Snow Model Performance............................................................................................... 94 3.2 Snowpack Analysis ......................................................................................................... 96 3.3 Recharge Analyses .......................................................................................................... 98 3.4 Statistical Analyses ....................................................................................................... 100 4. Discussion ........................................................................................................................... 100 5. Conclusions ......................................................................................................................... 102 Acknowledgments .................................................................................................................. 103 APPENDIX ................................................................................................................................. 104 CHAPTER 5: A WARMER FUTURE: CLIMATE CHALLENGES FOR MICHIGAN’S SKI INDUSTRY ................................................................................................................................ 106 Abstract ................................................................................................................................... 106 1. Introduction ......................................................................................................................... 106 2. Methods .............................................................................................................................. 112 2.1 Study Region................................................................................................................. 112 2.2 Data Sources ................................................................................................................. 114 2.3 Analysis Methods.......................................................................................................... 115 3. Results and Discussion ....................................................................................................... 116 4. Conclusions ......................................................................................................................... 124 Acknowledgments .................................................................................................................. 124 REFERENCES ........................................................................................................................... 125 viii LIST OF TABLES Table 2.1. Regional snowpack and snowmelt year type statistics. ............................................... 24 Table 2.2. Regional annual basin yield statistics. ......................................................................... 27 Table 2.3. Regional recharge statistics. ........................................................................................ 30 Table 2.4. Regional statistical analysis results.............................................................................. 31 Table A.3.1. Winter type statistics. ............................................................................................... 80 Table A.3.2. Statistical tests results. ............................................................................................. 81 Table 4.1. Regional winter type statistics for selected hydrologic model outputs. ...................... 97 Table 4.2. P-values for statistical tests. ....................................................................................... 100 Table 5.1. Winter type statistics. ................................................................................................. 118 ix LIST OF FIGURES Figure 2.1. Study region maps. ..................................................................................................... 11 Figure 2.2. HUC-8 basin multimetric scores. ............................................................................... 21 Figure 2.3. Median regional SWE in different year types. ........................................................... 23 Figure 2.4. Regional year type total melt, melt events and complete melts. ................................ 25 Figure 2.5. Median regional basin yield in different year types. .................................................. 26 Figure 2.6. Regional monthly net recharge and cumulative net recharge in different year types. 29 Figure A.2.1. Normalized seasonal temperature across the domain. ............................................ 40 Figure A.2.2. Winter PET for study region. ................................................................................. 41 Figure 3.1. Study region maps. ..................................................................................................... 47 Figure 3.2. Regional wam winters and metric scores per decade. ................................................ 55 Figure 3.3. Metric score time series. ............................................................................................. 56 Figure 3.4. Regional bare ground days maps. ............................................................................... 58 Figure 3.5. Regional winter type snow depths. ............................................................................. 60 Figure 3.6. Regional winter type basin yield. ............................................................................... 62 Figure 3.7. Regional winter type cumulative basin yield. ............................................................ 63 Figure 3.8. Regional spring vs annual discharge. ......................................................................... 64 Figure A.3.1. Additional study region map. ................................................................................. 73 Figure A.3.2. Regional winter type precipitation. ........................................................................ 74 Figure A.3.3. Regional winter temperatures. ................................................................................ 75 Figure A.3.4. Methodology flowchart. ......................................................................................... 76 Figure A.3.5. Regional metric score slopes. ................................................................................. 77 x Figure A.3.6. Regional winter vs annual discharge. ..................................................................... 78 Figure A.3.7. Regional basin yield difference. ............................................................................. 79 Figure 4.1. Study region................................................................................................................ 86 Figure 4.2. Selected inputs for LHM. ........................................................................................... 91 Figure 4.3. LHM snowpack model compared to SNODAS and GHCN data............................... 95 Figure 4.4. HUC-10 SWE totals for each complete water year in the study period. .................... 96 Figure 4.5. Mean monthly SWE in warm and cool winters.......................................................... 97 Figure 4.6. Mean simulated recharge by winter type.................................................................... 98 Figure 4.7. Cumulative monthly warm and cool winter recharge. ............................................... 99 Figure A.4.1. Winter type classification flowchart. .................................................................... 105 Figure 5.1. Study region map. ..................................................................................................... 112 Figure 5.2. Winter type daily depth. ........................................................................................... 117 Figure 5.3. Winter type depth and snowmaking suitability. ....................................................... 120 Figure 5.4. Snowmaking water use and hours. ........................................................................... 121 xi KEY TO ABBREVIATIONS CDV Center of Discharge Volume DEM Digital Elevation Model DJF December-January-February DOWY Day of Water Year ET Evapotranspiration GHCN Global Historical Climatology Network gSSURGO Gridded Soil Survey Geographic Database HUC Hydrologic Unit Code KS Kolmogorov-Smirnov LAI Leaf Area Index LHM Landscape Hydrology Model LP Lower Peninsula LULC Land Use-Land Cover MAM March-April-May MDEQ Michigan Department of Environmental Quality MODIS Moderate Resolution Imaging Spectroradiometer NED National Elevations Dataset NHD National Hydrography Dataset NLCD National Land Cover Dataset NLDAS North American Land Data Assimilation System NLP Northern Lower Peninsula xii NOAA National Oceanic and Atmospheric Administration NSAA National Ski Areas Association NWI National Wetlands Inventory PET Potential Evapotranspiration SLP Southern Lower Peninsula SNODAS Snow Data Assimilation model S/P Snow/Precipitation SWE Snow Water Equivalence Q Stream Discharge QA/QC Quality Assurance/Quality Control UEB Utah Energy Balance model UP Upper Peninsula USDA US Department of Agriculture USGS US Geological Survey xiii CHAPTER 1: INTRODUCTION As the climate continues to warm (IPCC, 2022) it is changing the timing and amount of snowfall in alpine and mid-to-high latitude regions across the globe. Many of these areas rely on the water produced from these reduced snowpacks, and the changes brought about by the warming are not entirely obvious. However, the more these snowpacks’ annual patterns deviate from the norm, the more significant the possible effects on winter and spring hydrology as streamflow and groundwater patterns shift in response (Hyndman, 2014). Much of the research into changing snowmelt dynamics from climate warming are focused on the mountainous regions of the planet where the snowpack constitutes a larger portion of the annual water budget, but many non-alpine, temperate high latitude regions of the Earth such as the Great Lakes Basin receive substantial snowfall amounts every year. With climate change continuing to progress unabated, snowpack and snowmelt will become critical for subsequent hydrologic changes in these regions. The Great Lakes region receives upwards of 300 cm/year of snowfall in the northern parts, with most of that snow falling from October through May (Ford et al., 2021). Overall annual precipitation amounts average 82.6 cm, with 49.5 cm of that precipitation falling during the October-May snow season. In addition to the cold winter temperatures in these northern states, much of this snow in the Great Lakes is driven by the “Lake Effect” phenomena, where the large, open bodies of water of the lakes experience increased evaporation, creating moisture- laden clouds which then drop much of this moisture as rain or snow shortly after reaching land. Typically, when these heavy snowfall amounts occur, the snowpack accumulates on the ground throughout the winter and early spring before melting late in the spring, generating large 1 amounts of meltwater that recharges shallow groundwater and causes high spring streamflow. However, with the winter temperatures in these regions becoming warmer, this seasonally persistent snowpack is becoming more episodic, with thinner snowpacks that melt (in some cases entirely) throughout the season, causing “pulses” of meltwater through the winter and spring rather than one large melt at the end of the season. This change to the pattern of melting could significantly alter the winter and spring hydrology of these regions, changing how much of that melt generates streamflow vs. how much melt infiltrates the subsurface and recharges groundwater. The challenges in studying snowpack and snowmelt dynamics stem from the inherent heterogeneity of the snowpack in both space and time. Spatially, snowpack thicknesses can vary substantially at very small spatial scales. The combined effects of the landscape topography, wind movement and vegetation cause significant changes in thickness to occur within small distances. This makes traditional methods of snow data collection from individual weather stations unrepresentative of the general snowpack thickness. This is further complicated by the evolution of the snowpack through time. Oftentimes the snowpack can melt with that meltwater refreezing within the snowpack, causing layering of the snowpack with dense, icy layers near the bottom and fresher, less dense snow layers near the top. As the snowpack matures, the snow becomes denser, decreasing the snow depth while losing little water content. This is a large part of why snow depth is a less reliable measurement of the water content of the snow, and for hydrologic research the snow-water equivalence (SWE) is a better variable. The SWE is the amount of water content held by the snow, but unfortunately it is a difficult snow parameter to measure in the field. Thus, most snow hydrology researchers rely on modeling, both statistical and simulation, to better grasp snowmelt hydrology. 2 Several methods for modeling snowpack progression and snowmelt generation exist. Perhaps the easiest modeling approach is to aggregate the point-scale observational data from weather stations across large spatial and temporal domains to detect statistical trends in the snow. Snow stations can be grouped by contiguous, regional shapes such as drainage basins, and then average depths can be extracted over daily to annual timesteps. The limitation of this approach is the spatial resolution because regions grouping stations together ideally are large enough to capture multiple stations, but the larger this area the less accurate it will be in capturing the smaller-scale variations of the snowpack and may not be a good representation of the snowpack as whole. To capture those local snowpack variations, the snowpack needs to be modeled in a simulation. With simulation modeling the temporal and spatial resolutions are limited only by time and processing power of the hardware, allowing examination of small-scale snowpack heterogeneities that isn’t possible with empirical data. However, a simulation model is only as accurate as each of its parts, and the model requires a complex interaction of inputs relating to precipitation, energy and heat fluxes, geology, and hydrology. If any of these model pieces are inaccurate, the entire simulation’s accuracy is suspect. When properly calibrated, though, the snowpack and the melt produced from it can be examined at finer scales, allowing for more thorough investigations into the snow’s contribution to the local water budget. Assessing these changes to local and regional hydrology in non-alpine settings resulting from changes in snowmelt was the primary goal of this dissertation research, addressing current gaps in the scientific literature while thinking about the real-world implications of these changes. While the research in this dissertation focused on the Great Lakes basin, spatial and temporal scales varied for each study, with some focusing just on the state of Michigan at 1–2-decade 3 temporal periods and others looking at the Great Lakes and eastern United States over half a century. The methodology of each study also evolved over the course of the dissertation, starting with the observational weather station data and modeling outputs from federal scientific agencies to complex simulation of the landscape hydrology using modeling products locally developed. Despite the differences in scales and approach, the individual studies of this dissertation each address the changing winter climate in the Great Lakes and what that means for surface and groundwaters in the region. Chapters 2 and 3 of this dissertation focus on the processes of changing snowmelt in the Great Lakes, relying on empirical data for establishing that winters are in fact becoming warmer with time and how the snow and hydrology are responding to this warming. Chapter 2 uses a combination of observational point-scale data and modeled outputs of precipitation and the snowpack from 2003-2017 to identify winters as warmer or cooler than the norm during this period and observe the differences in melt and hydrology during these different winter types. Chapter 3 then takes the methodology established in the previous chapter and implements it across much of the eastern United States from 1960-2019, defining just how much warmer winters are getting and the extent of the hydrologic changes. Chapter 2 is published in the Journal of Hydrology (Ford et al. 2020) and Chapter 3 is published in Science of the Total Environment (Ford et al. 2021). Both chapters have been reprinted here. Using these findings, Chapter 4 simulates the changes to snow and hydrology in the Landscape Hydrology Model (LHM) across the state of Michigan from 2000-2019. The model captures the observed changes to melt and streamflow during this period while filling in the critical missing component of groundwater recharge, which is difficult to examine due to the lack of observational data regarding groundwater levels. The results of the simulation then 4 inform the final chapter of this dissertation, Chapter 5, which examines the potential economic impacts of these changes in the region with specific focus on how changes to the snowpack in warmer winters may affect the industries reliant upon that snow such as ski slopes and other winter tourism. In total, this research fills in scientific knowledge gaps regarding changes to snowmelt dynamics in non-alpine temperate settings in response to warming winter temperatures. In addition to establishing the shifting patterns of snowpacks and melt-driven hydrology, these studies examine the often-overlooked resulting groundwater changes. These insights into the changing snowmelt hydrology will allow better forecasting of future shifts as the climate continues to warm. Such information will be critical to resource managers and ecological stewards in the future as humanity grapples with the changing environment around them and determines best practices moving forward in an uncertain future. 5 CHAPTER 2: EFFECTS OF SHIFTING SNOWMELT REGIMES ON THE HYDROLOGY OF NON- ALPINE TEMPERATE LANDSCAPES Abstract Across the globe, increasing winter temperatures resulting from climate change are affecting how and when snow melts. Despite the significance of snowfall to the annual water budget in mid-latitude, non-alpine regions, there has been limited research into shifting snowmelt dynamics in these areas. This study examines snowmelt changes resulting from warmer winters in the Great Lakes region of North America. Using multiple metrics calculated from station temperature data, recent years (2003-2017) were categorized as “warm” or “cool”, then snowmelt and regional hydrologic patterns within those year types were examined. Systematic differences were observed, including warmer years having less and earlier snowmelt than cooler years. Those changes have, in turn, led to both lower and earlier spring peak flows in streams and decreased net groundwater recharge in the northern regions. Additionally, we show that differences between warm and cool year types become more significant along a north-south gradient; differences between warm and cool years are more pronounced in northern regions with regard to streamflow, net recharge, melt amount and timing. These results provide a framework to examine linked changes between snowmelt and hydrology across high-latitude temperate regions of the world. 1. Introduction High-latitude, temperate, low elevation regions around the globe depend on seasonal snowfall for a significant portion of their annual water budgets. While cryosphere research has generally focused on alpine and arctic regions, changing climate conditions are already affecting snow hydrology across vast but relatively understudied regions (Burakowski et al., 2008; 6 Hodgkins and Dudley, 2006a; Huntington et al., 2004; Javed et al., 2019; Suriano et al., 2019). Indeed, the comparatively low-relief topography characteristic of the non-alpine, non-arctic cold regions of the world make large areas particularly sensitive to the effects of shifting patterns of snow accumulation and melt (Clark et al., 2011). In particular, these regions are moving from relatively persistent seasonal snowpack to a thinner snowpack and increased bare ground days (Suriano and Leathers, 2017). The hydrologic effects of this regime shift are not well understood, but will reshape economies, ecosystems, water resources, and the needs of built infrastructure during the coming decades (Chin et al., 2018; Suriano et al., 2019; Zou et al., 2018). Snowpack and melt research in recent decades has primarily focused on the mountainous regions of the western United States, Europe, and Asia. Many studies have examined the potential effects of global climate change on alpine snowpacks and subsequently on mountain streamflow. Studies of historical trends in snow accumulation and melt in the western United States have shown declining snowpack thicknesses, earlier melting, and more winter/spring precipitation occurring as rain rather than snow over recent decades (Clow, 2010; Hamlet et al., 2005; Jefferson et al., 2008; Mote, 2003; Stewart et al., 2004). These trends are not limited to the alpine areas of the country. Several regional studies have found similar historical snow trends in the central and eastern U.S. (Burakowski et al., 2008; Hodgkins and Dudley, 2006a; Huntington et al., 2004), including reduced snow cover, lower snow to total precipitation ratios (S/P) and more melting in winter months rather than the traditional late Spring melt (Dyer and Mote, 2006; Feng and Hu, 2007; McCabe and Wolock, 2009; Suriano and Leathers, 2017; Suriano et al., 2019). Future climate scenarios continue to project decreased snow with earlier melting (Adam et al., 2009; 7 Barnett et al., 2005; Boyer et al., 2010; Campbell et al., 2011; Chin et al., 2018; Demaria et al., 2016; Hayhoe et al., 2007; Hayhoe et al., 2010; Mortsch et al., 2000; Peacock, 2012). The hydrologic effects of these melt changes have been observed in alpine settings of the western U.S., including lower peak flows that occur earlier in the season, and earlier overall streamflow regimes including earlier center of streamflow volume arrival times (Clow, 2010; Hidalgo et al., 2009; Jefferson et al., 2008; Stewart et al., 2009). Similar observations of spring streamflow have been made in the eastern half of the country (including the Great Lakes), but to a lesser degree (Hodgkins and Dudley, 2006b; Hodgkins et al., 2007; Huntington et al., 2004; Javed et al., 2019; Johnson and Stefan, 2006; Johnston and Shmagin, 2008). As the warming associated with global change progresses, the expectation is that streamflow will continue to shift earlier in the year across the globe (Arora and Boer, 2001; Barnett et al., 2005; Berghuijs et al., 2014; Boyer et al., 2010; Brubaker and Rango, 1996; Byun and Hamlet, 2018; Campbell et al., 2011; Champagne et al., 2020; Cherkauer and Sinha, 2010; Hayhoe et al., 2007; Mortsch et al., 2000; Tu, 2009). As such hydrologic changes become more pronounced and widespread, it will be critical for research to be conducted on the snowmelt hydrology of lesser studied regions of the planet, including the often-overlooked groundwater component of the hydrologic cycle. Here we focus on the Laurentian Great Lakes region of North America (hereafter “Great Lakes”), and in particular the state of Michigan, USA. Situated in the midst of four Great Lakes, Michigan is geographically unique, and serves as a mesocosm for high-latitude temperate climate regions globally. The range of snowfall across the state, both due to its latitudinal extent and the “Lake Effect” phenomenon, encompasses nearly the entire range for the continental US; some areas within the state receive less than 100 mm SWE of snowfall, while the northern regions typically receive just over half a meter in a year (Figure 2.1a). Michigan also has 8 remarkable hydrogeologic diversity for its area, spanning 12 of the 20 US Hydrologic Land Regions (Wolock et al., 2004), with soils spanning the range from heavy clay to coarse sand. Despite the importance of snow to Great Lakes annual water budgets, and the potential negative consequences of changes to the typical snowmelt regime, research into snowmelt patterns across this region has been limited. Within the Great Lakes region, relatively few studies have examined historical snowmelt trends. Hodgkins et al. (2007) found that from 1955-2004 snowmelt occurred earlier in the year with decreasing S/P ratios. Most climate model simulations of the Great Lakes region show increases in winter and spring precipitation, however warmer temperatures will lead to more rain and less snow (Hayhoe et al., 2010). Historical data from the last century in this region shows that as climate warmed, annual precipitation amounts have increased, leading to increased streamflow (Hodgkins et al., 2007). The relation of this increased precipitation to streamflow from historical data remains unclear, as early spring precipitation typically decreased over the last 50 years, but total runoff (sum of surface and groundwater components) increased, indicating earlier melting and more precipitation falling as rain. The decrease in amount of runoff from melt corresponded with flow peaks earlier in the season, decreased peak streamflows, and earlier peak Great Lakes levels (Argyilan and Forman, 2003; Johnson and Stefan, 2006; Novotony and Stefan, 2007). Those studies analyzed changes to the hydrologic cycle focused on a specific drainage basin (or similar scales). Where studies have taken a broader spatial perspective, they have generally focused on analyzing just one component of the hydrologic cycle. Here, we seek to holistically examine how changing snowmelt dynamics are affecting the broader hydrologic cycle over a large non-alpine region. To best understand the potential snowmelt hydrology changes, taking into account differences in climate, land use, geology and hydrology across Michigan, we quantify changes to 9 seasonal snowmelt and the influence it has on groundwater recharge and streamflow. Anecdotally, Michigan’s seasonal snowpack in Michigan in recent years appears to have deviated from the historically “typical” season-long persistence, instead melting periodically throughout the winter months. We hypothesize that: 1) warmer winters have led to less snowmelt, and 2) reduced snowpack persistence, typified by earlier melt and more bare ground days throughout the winter, leading to 3) earlier and lower spring peak flow in streams, and 4) increased recharge of shallow groundwater. This paper examines these four linked hypotheses by first categorizing recent years as “warm” or “cool” based on a multimetric analysis of different winter temperature parameters across Michigan’s substantial north-south climate gradient. We leverage 14 years of gridded, assimilated model data of snowpack and melt to quantify differences across this gradient between warm and cool years. Finally, we correlate these changes with variations in drainage basin hydrology across year types. 2. Methods 2.1 Study Region Although the high-latitude Midwestern United States, particularly the Great Lakes region, is neither arctic nor alpine it still receives a significant portion of its water budget from seasonal snowfall. Some of the highest annual snowfall totals in the eastern half of the United States occur in this region. In the Upper Peninsula (UP) of Michigan, annual snow water equivalence totals can exceed half a meter (Figure 2.1a, and Andresen, 2012). 10 Figure 2.1 Study region maps. Study region and data source maps showing Michigan in the inset map, along with: a) annual average snow water equivalent from 2003 – 2016 computed from SNODAS reanalyses overlain by the three generalized regional polygons used in snowmelt analysis; b) Global Historical Climatology Network (GHCN) station locations, along with the HUC-8 basins; c) available USGS gauge locations in Michigan and their drainage basins, colored by three study-defined sub-regions; d) the Land Use-Land Cover across the state from the National Land Cover Dataset (USGS, 2011); e) the heavily glaciated surficial deposits (EGLE, 2019); and f) the saturated conductivity of the top soil layer from gSSURGO (Soil Survey Staff, 2020). Impacts of warming winter temperatures due to global climate change on Michigan’s winter hydrology are unclear. Climate projections for Michigan show an increase in annual precipitation, but it is not clear how much of that precipitation will occur as snow (Hayhoe et al., 2010); even in areas where annual snowpack thickness has been increasing, the number of days 11 with snow covered ground appears to be decreasing (Andresen, 2012). Fewer days with snow on the ground would indicate a shift in the melt regime of the snowpack, in turn altering recharge and spring streamflow amounts. Michigan’s geology changes along a north-south gradient due to numerous glacial advances and retreats during the last ice age (Figure A.2.1; Farrand and Bell, 1992; Groundwater Inventory and Mapping Project, 2003). The shallow sediments and aquifers in the LP are dominantly sands to sandy loams, which are mostly underlain by sedimentary bedrock. In contrast, the UP has shallow soils underlain by less permeable Precambrian igneous formations. As a result, it has distinctly different hydrology that is more dominated by surface flow than subsurface flows. Michigan’s also has a strong land use gradient from north to south. The southern portion is predominantly agricultural, growing the same variety of staple crops seen throughout the rest of the U.S. Midwest: corn, soybeans, alfalfa and hay crops are some of the primary agricultural outputs of the region (Hamlin et al., 2020; Homer et al., 2004; 2015; Figure A.2.1). Most of the state’s major urban centers are also in the south, creating a landscape of urbanized areas surrounded by intensively managed farmlands. By contrast the state’s northern Lower Peninsula (NLP) and Upper Peninsula (UP) are dominated by mixed and coniferous forests with less urban and agricultural land than the southern portion. There are distinct precipitation patterns across Michigan, largely due to the “Lake Effect” phenomenon (Figure 2.1). The LP’s climate is strongly influenced by winds that come from the north- to the south-west across Lake Michigan (Andresen, 2012). As the air moves over the large lake, water vapor content increases due to evaporation. Shortly after these air packages reach land, this additional water vapor commonly condenses, producing substantially more 12 precipitation on the western portion and northern tip of the peninsula. Average annual snow water equivalent (SWE) totals across the LP range from < 200 mm in the southeastern LP to >400 mm in the northern part. The UP’s precipitation patterns are similarly affected by Lake Superior. The northern parts of the UP can receive very heavy snowfall, with an average annual SWE of up to ~ 0.5 meter. To best understand potential snowmelt hydrology changes, taking into account differences in climate, land use, geology and hydrology, the state was divided into three regions (Figure 2.1). The UP, with its high wetland and forested fraction and distinct hydrogeology was assigned as one region. The LP was split in two (Northern Lower Peninsula, NLP, and Southern Lower Peninsula, SLP), with the dividing line of approximately 43.8° N latitude chosen because it best divides the LP based on geologic, land use, and climatic differences. Using these regional delineations, snowmelt patterns and hydrologic responses were analyzed between the 2004 to 2017 water years. 2.2 Data Sources We combined observations and model-data reanalysis from several sources. Daily temperature and precipitation data from weather stations were extracted from the Global Historical Climatology Network (GHCN) (Menne et al. 2012). Daily snowpack SWE and melt values were derived from the National Oceanic and Atmospheric Administration’s (NOAA) Snow Data Assimilation (SNODAS) model (NOHRSC, 2004). Stream gauge daily average flow data came from the U.S. Geological Survey’s (USGS) stream gauge network (USGS, 2018). Daily precipitation and average temperature values were extracted from the PRISM reanalysis data (PRISM Climate Group, 2018). All datasets were analyzed from October 2003 to June 2017. 13 All GHCN and USGS gauges within the three study subregions were included, provided that 95% of the total daily observations across “snow seasons” (defined here as October 1 through May 31 of the following year) were present. With this constraint, maximum and minimum daily air temperature were available from 240 weather stations. Mean air temperature for each day was calculated as the arithmetic mean of the provided maximum and minimum daily temperatures. Daily streamflow data were available from 123 USGS stream gauging stations. The output from NOAA’s SNODAS model was one of the primary input sources for this study. SNODAS is a data assimilation model that is calibrated to snow and climate observation data (Barrett, 2003); it performs well, particularly in areas with relatively low relief (Clow et al., 2012; Hedrick et al., 2015). SNODAS uses downscaled outputs from numerical weather prediction models in conjunction with empirical meteorological data from airborne and ground- based weather stations along with satellite data to model snow across the continental United States. The model outputs daily 1-kilometer grids of snowpack thickness and temperature, along with other simulated snowpack properties such as daily melt, across the conterminous United States. The model output first became available in early 2003, which is thus the beginning date for our analysis period. We used the outputs of modeled snowpack SWE and melt from the base of the snowpack. The non-snow precipitation data came from the PRISM model (PRISM Climate Group, 2018). This model outputs daily 4-kilometer gridded interpolations of total precipitation, minimum and maximum temperature. Data was downloaded through FTP using the R package ‘prism’ (Hart and Bell, 2015). Since the model doesn’t differentiate between rain and snow, we assumed that all precipitation fell as snow below a threshold of 1.5° C, and as rain above this 14 threshold. This threshold is near the upper bound of the transition temperature across non-alpine North America found in the study by Jennings et al. (2018). 2.3 Analysis Methods Daily input data were aggregated spatially and temporally to classify years as “warm” or “cool”, to analyze changes in snowpack and streamflow, and compute basin-wide water budgets. First, we developed a multimetric analysis of the GHCN temperature data to classify years into “warm” and “cool” relative to the mean for the 14-year period. Then, using these yearly classifications, SNODAS snow and USGS streamflow data were analyzed within each year type to evaluate differences in melt and streamflow amounts and timing. The melt estimates from SNODAS were used in conjunction with the streamflow data to evaluate hydrologic effects of melt changes. Statistics of seasonal flow timing and amount were extracted from the daily gauge data across the basin. These streamflow data were then compared to SNODAS output and precipitation data from PRISM to examine seasonal net recharge associated with melt. All of the data analysis techniques were performed in R. 2.3.1 Spatial Aggregation Several spatial aggregations were employed during the analyses. The GHCN station data were aggregated over HUC-8 basins to provide complete spatial coverage of the state and allow temperature variation among basins. The SNODAS spatial data were imported as raster files into the statistical software R (R Core Team, 2019) and then mean daily melt/SWE values were extracted across stream basin polygons, which were calculated for each USGS gauge station using 30 meter (1 arc second) National Elevations Dataset DEM data. The DEM-based stream basin polygons were used for this analysis to relate streamflow to the snowmelt and precipitation inputs over each basin. These spatial aggregates across basins were then viewed through the lens 15 of the three broader regions of the state (SLP, NLP and UP), which are used to analyze melt, streamflow and net recharge. 2.3.2 Multimetric “Warm” and “Cool” Year Classification We developed a multimetric classification of winter/spring air temperatures as “warm” or “cool” to go beyond simple seasonal average temperatures typically used. We first aggregated the station data spatially, computing daily average air temperatures within HUC-8 basins (see Figure 2.1). Hydrologic Unit Code (HUC) basins are drainage basins delineated as part of the “Hydrologic Unit Maps” project of the USGS, ranging from the largest two digit HUC’s (HUC- 2’s) to the smallest 12 digit basins (HUC-12’s) (Seaber et al., 1987; USGS, 2014). HUC-8 basins were used to group stations because they offer complete coverage of the state without any overlap, generally include one or more temperature stations, and are sufficiently small to not obscure smaller-scale variations in temperatures. Within each basin, we then computed water- year values for six air temperature metrics: 1) average from October-May; 2) minimum for the same period; 3) winter (Dec-Feb) average and 4) winter minimum; 5) spring (Mar-May) average, and 6) spring minimum. These metrics were chosen because temperature changes have been found to be the primary driver for changing snowmelt dynamics (Boyer et al., 2010; Cline, 1997; Hamlet et al., 2005; Hodgkins and Dudley, 2006a, 2006b; Hodgkins et al., 2007; Jefferson et al., 2008; Johnson and Stefan, 2006; McCabe and Wolock, 2010; Mote, 2003; Stewart et al., 2004). We sought to include metrics that would capture whole-season temperatures, equally and separately weighting both winter and spring portions. Minimum temperatures were included because if the nighttime air temperatures stay above the melting temperature the snowpack continues to melt rather than refreeze; minimum temperatures below freezing would delay 16 melting the following day as warmer daytime temperatures must reheat the snowpack until an isothermal snowpack temperature gradient can be reached to generate melt (Cline, 1997). “Warm” and “cool” year classifications were then determined from multimetric z-scores. Each metric’s z-score was calculated by subtracting the yearly mean and dividing by the 14-year standard deviation within HUC-8 basins. The arithmetic mean of the six resultant single metric z-scores (each centered around 0) was then calculated to provide a single annual multimetric score for each HUC-8 (Figure A.2.1). An overall annual state metric score was then calculated as the mean across all HUC-8 basins. More positive values for this overall annual score indicate warmer winters than the 14-year norm, and negative values thus indicate cooler. Water years with multimetric scores less than -0.5 were classified as “cool” and greater than 0.5 were classified as “warm”. Metric scores between -0.5 and +0.5 were deemed “normal” years. For an individual metric, the 0.5 threshold represents 0.5 standard deviations away from the mean; for the multimetric score the values have a similar meaning though they are not precisely defined as standard deviations above or below the norm. These classifications were then used as the basis for the remainder of the analysis in this study. 2.3.3 Regional-Scale Snowpack and Streamflow Analyses After classifying each year in the study period as warm, cool or normal the SNODAS data were analyzed within each of our three regions (UP, NLP, and SLP) for changes to melt amount and timing. Daily mean SWE and mean melt output were averaged across grid cells within each of the three regional polygons. Then, annual values for the following statistics were calculated for each region: peak melt/SWE amount, peak melt/SWE timing, number of bare ground days, annual melt amount, 50th quantile of seasonal melt volume (SM50), number of melt events (defined here as periods of consecutive days with melt generated from the snowpack), 17 melt event length and amount, and number of complete melt events (defined as melt events that ended with no remaining snowpack). We then computed arithmetic means of each of these statistics within warm and cool year types. We also computed warm and cool year average daily time-series SWE curves within each region. These were calculated as the average day-of-water-year SWE within each regional polygon for each year type. This provided an informative visualization of the general snowpack progression through the season, and illustrates both regional and year type differences within the snowpack. We examined the distribution of the data using several different methods. Most figures of amounts and timing across station/gage data were plotted as violin plots; these plots display y- axis values similar to a standard boxplot, but also display the density of data around a given y- axis value as a vertically-mirrored kernel density estimate. This shows the distribution of data along with the quantiles. Finally, we compared year type distributions using two statistical tests described below. Stream gauge data were analyzed in a similar manner as the SNODAS output. Since these gauges provide point data that are associated with individual drainage basins for each gauge, there was no need for any spatial aggregation. For each gauge site, statistics were calculated for: winter season peak flow amount and timing, annual flow quantiles, and basin yields. Basin yield is defined as the daily streamflow divided by gauge basin area. Values were then averaged across all gauge basins within each of the three regions. Similar to the snowpack time-series SWE curves, we computed a mean basin yield hydrograph across regions within warm and cool years. From these, we then computed the center of volume (CV), defined as the date of arrival of 50% of the flow between October and May. 18 2.3.4 Winter and Spring Net Recharge Analysis Within each stream gauge basin we computed the daily overall basin water balance. This water balance allowed us to estimate net winter and spring recharge (i.e. the net change in storage). Assuming there is minimal evapotranspiration between December and April (Kirchner and Allen, 2020; Figure A.2.2), and that streamflow is the only significant basin outlet (i.e., pumping or losses to deeper geologic units are minimal) the water balance for each gauge basin can be written as ∆𝑆 = 𝑃 + 𝑀 − 𝑄 (1) where ΔS is the change in groundwater storage (i.e., net recharge, mm/d), P is rain (mm/d), M is snowmelt (mm/d) and, Q is stream discharge expressed as daily basin yield (mm/d). Rainfall was computed from PRISM precipitation as described above, and daily snowpack basal melt data were extracted from SNODAS to quantify total available liquid water (𝑃 + 𝑀). Daily averages across grid cells for PRISM and SNODAS were calculated within each stream gauge basin. To reduce noise, these daily values were then summed monthly, and medians of the monthly sums was calculated across regions and year types. 2.3.5 Dataset Distribution Tests Two statistical tests were used to evaluate if data and calculated values from warm years were statistically different from those in cool years. First, annual basin values such as total seasonal melt, the timing of peak SWE, total rain and peak basin yield were calculated. These annual datasets were classified by their year type, and then warm and cool year distributions were evaluated for normality using the Shapiro-Wilk Normality Test (Royston, 1995). These year type data were then compared using both the two-sample Mann-Whitney (Wilcoxon) ranked sum test (Hollander and Wolfe, 1973) and the two-sample Kolmogorov-Smirnov (KS) test 19 (Conover, 1971). The p-values for these two tests were evaluated on a 95% confidence threshold to evaluate if the two year-type datasets were statistically different. 3. Results 3.1 Warm and Cool Year Classification Four years were classified as “warm” (2004, 2010, 2012, 2017) and five years as “cool” (2008, 2011, 2013, 2014, 2015) using the multimetric analysis. The spatial distribution of metric scores is generally uniform (Figure 2.2) across the study area. In particular, spatial homogeneity increases as metric scores deviate from 0. A visual examination of the metric score distribution does not reveal any temporal trends, with the coldest year (2014) and warmest year (2012) occurring within a short timeframe. Nor does the relatively short study period provide a long enough timeframe to robustly assess trends. The influence of the individual metrics on the overall metric score was examined by comparing time series of the individual metric scores averaged across all HUC-8 basins in the study area (Figure A.2.1). While there was deviation from the overall metric score for some metrics, rarely did individual metrics classify years differently from the overall metric. The average spring temperature and minimum spring temperature metrics classified the year as cooler than the overall metric two and three times respectively. The winter average and minimum temperature metrics only classified a year as cooler than the overall metric once. 20 Figure 2.2. HUC-8 basin multimetric scores. Annual maps of normalized multimetric scores for each HUC-8 basin from water year 2004 to 2017. White HUC-8’s are basins without data from a GHCN station. Panel titles are colored to denote year type based on statewide multimetric score: red for warm years and blue for cool. 3.2 Precipitation Analysis Using the K-S test on the PRISM data, we examined if the hydrologic differences between year types were driven by changes in total precipitation (both solid and liquid) amounts rather than changes to snowmelt regimes. The K-S test was run on total precipitation as well as each of the precipitation types individually comparing warm year precipitation totals to cool year totals. Gridded total precipitation across the study region showed significant differences between 21 cold and warm winters (P<0.001), with warm winters having an average of 356 mm and cool years 419 mm. While both snow and rain distributions between the two year types have P-values below 1% (P < 0.001 and P = 0.008, respectively). 3.3 Snowpack Analysis Snow in all regions generally starts to accumulate in late November (around the 50th day of the water year, Figure 2.3), with the UP receiving its first snow on average 21 days before the Southern Lower Peninsula (SLP). Cooler winters have consistently higher SWE regardless of the region or time of year, and the persistence of the snowpack is significantly longer. In both year types, snow persists significantly longer in the UP than in either the Northern or Southern LP (on average lasting 17 and 20 days longer than snow in the Northern LP in warm and cool years respectively). Across regions and year types average peak SWE occurs in a 6-week window, ranging from January 28-March 9. Snow amounts increase moving northward for all year types, with the UP’s maximum SWE on average about 100 mm higher than the SLP maximum SWE (Figure 2.3, Table 2.1). Note that snowpack SWE variability is highest between days 160 and 200 (mid-March to late April). 22 Figure 2.3. Median regional SWE in different year types. Median daily snow-water equivalence (SWE) for each year type for a) the UP, b) the Northern LP and c) the Southern LP. Cool years (blue) show higher SWE and persist later into the spring than warm years (red) in all regions. Time is plotted as day of water year starting on October 1. Upper and lower bounds of the shaded regions represent the maximum and minimum, respectively. These regional differences that are independent of year type also appear in the total seasonal melt within HUC-8 basins, with those in the UP experiencing nearly double the amount of melt as the SLP in both year types (Figure 2.4). These regional differences are also reflected in winter bare ground days (Table 2.1). In all year types, the UP has fewer bare ground days, correlating with fewer complete melt events. The SLP has the most bare ground days regardless of year type, while the NLP falls between the two. The snowpack in the UP melts in fewer melt events than the more southern regions, and relative to the SLP, most of those melts do not continue to complete melt of the snowpack. The NLP shows differences and similarities to the other two regions, having similar numbers of melt events per season as the UP, but more of those melt events continue to completion, similar to what’s seen in the SLP. 23 Region Southern LP Northern LP UP Year Type Cool Warm Cool Warm Cool Warm Peak SWE (mm) 72 35 99 74 167 123 Peak SWE Day of Water 132 119 132 147 159 153 Year Bare Ground Days 128 157 90 135 60 100 Peak Melt (mm) 38 23 38 40 40 29 Peak Melt Day of Water 121 140 158 167 217 174 Year 50% Melt Day of Water 157 140 177 163 210 180 Year Total Melt (mm) 304 215 420 296 581 420 Melt Events 9 8 7 8 5 7 Melts to Completion 5 4 2 5 2 4 Melt Event Amount (mm) 24 20 15 10 26 13 Melt Event Length (days) 9 7 10 7 24 10 Table 2.1. Regional snowpack and snowmelt year type statistics. Regional snowpack and snowmelt statistics averaged within year types. Bare ground days are calculated from October 1 to May 31. More melt is produced (and thus more snowfall occurred) in cool years regardless of region, with cool years producing around 100 mm more melt in all regions (Table 2.1). The Mann-Whitney ranked sum test on annual basin melt totals grouped by year type showed a significant difference between melt in the two different year types with a p < 0.0001 (Table 2.4). Cool years have more total melt events but fewer that go to completion where no snow is left on the ground (Figure 2.4, Table 2.1). In warm years in all regions the timing of peak melt closely correlates with the timing of 50% of seasonal melt, but in cool years, peak melt occurs several weeks before this date. Regionally these year type differences are greater in magnitude in the UP and SLP. These two northern and southern regions show the most significant differences in snowfall and snowmelt between year types, while the year type differences in the NLP are less significant (Table 2.4). The distribution of snow-related values between warm and cool years is significantly different with a 95% confidence interval in all regions except the NLP, which had p-values 24 below this confidence interval for peak melt amount and timing, maximum SWE timing and the timing of 50% of bulk seasonal melt (Table 2.4). In the UP, only two values were not statistically different according to the ranked sum test: total seasonal melt (p = 0.15) and season length (p = 0.077). In the SLP only one metric scored below the confidence interval threshold: peak melt time (p = 0.13). When comparing state-wide distributions between year types all values were significantly different for both tests. Figure 2.4. Regional year type total melt, melt events and complete melts. Violin plots of annual HUC-8 basin snowmelt statistics across regions (columns) and year types (blue is cool, red is warm). Rows are (top) total seasonal melt, (middle) number of melt events, and (bottom) ratio of complete melts to the total number of melt events. Here, the violin plot displayed the mirrored kernel density estimate along the y-axis, along with a horizontal bar for the median. There are large differences in total melt across regions and year types, with more melt events in warm years in the south. Additionally, within the SLP and NLP more melts in warm years continue to completion. While the SLP and UP snow datasets both show different distributions between year types, the way those datasets differ depends on the region. During the onset of winter in October- 25 December, the SLP snowpack shows similar patterns in both warm and cool years, but these differences increase as the season progresses (Figure 2.3). The median day of peak melt in the SLP is approximately 19 days earlier in cool years, but it is 43 days earlier in warm years in the UP (Table 2.1). It must be noted that this may be influenced by the small number of years in the study period. 3.4 Streamflow Analysis Figure 2.5. Median regional basin yield in different year types. Median daily basin yield (streamflow divided by basin area) across years and basins for the UP (a), Northern LP (b) and Southern LP (c). Days are plotted as days of the water year, beginning on October 1st. Vertical lines represent the center of volume day of water year for that year type during the October-May study period (denoted by the gray box). Similar to the snowmelt results, the differences between year types is largest in the SLP and UP (Figure 2.5, Table 2.4). The timing of the peak regional basin yield is earlier in all regions in warm years compared to cool years, ranging from one day earlier in the SLP to 36 26 days earlier in the UP. The center of volume timing is also earlier in warmer winters in all regions, most notably 30 days earlier than in the UP. The UP’s peak basin yield in cool years is twice that of peak yield in warm years (3.0 vs. 1.6 mm/day), but in the SLP, year type peak yield are similar with a difference of only 2.0 mm/day. Cool years produced more total seasonal basin yield in all regions, but the difference is only large in the UP. Southern Lower Northern Lower Region Upper Peninsula Peninsula Peninsula Year Type Cool Warm Cool Warm Cool Warm Total Basin Yield (mm) 200 199 231 222 194 161 Peak Basin Yield (mm) 2.0 2.2 1.6 1.9 3.0 1.6 Peak Basin Yield Day of 166 165 196 165 205 169 Water Year Center of Volume Day 164 147 140 136 178 148 of Water Year Coefficient of Variation 45 39 23 24 69 37 (%) Table 2.2. Regional annual basin yield statistics. Regional annual basin yield (discharge divided by basin area) statistics extracted from daily medians across year types and region. Warm years in the north have less total and peak basin yield, as well as earlier peak flow and center of volume date (day when 50% of that stream’s volume has occurred). Coefficient of variation is the ratio of the standard deviation of basin yield to the mean. Regional differences within year types are more noticeable in cool years. Peak basin yield in cool years increases northward, correlating with more melt ranging from 2.0 mm in the SLP to 2.0 mm in the UP. In warm years these regional differences are muted, with the peak basin yield increasing southwards in warm years, from 1.6 mm in the UP to 2.2 mm in the SLP (Table 2.2). The results from streamflow in the NLP don’t show as clearly defined a trend, with the NLP having the highest total basin yield amounts across regions in all year types, but the lowest peak basin yield amounts. 27 Once again, the lack of difference between year type datasets in the NLP is reflected in the results of the statistical tests run on the data distributions. Across the entire state, warm and cool year total seasonal basin yield, peak basin yield amount and timing, and the timing of 50% of seasonal cumulative flow (BY50) are all significantly different using a 95% confidence interval threshold for both the ranked sum test and the KS test (Table 2.4). Similarly, when examining the SLP and UP these basin yield datasets are all significantly different for both tests using this threshold. However, the NLP falls below the confidence interval in BY50 (ranked sum p = 0.32; KS p = 0.30), total seasonal basin yield (KS p = o.22) and peak basin yield timing (ranked sum p = 0.60). 3.5 Winter and Spring Recharge Analysis Net groundwater recharge (from Equation 1) monthly values across regions and year types show striking patterns (Figure 2.6). There were significant variations in net groundwater recharge estimates by region. Groundwater storage decreased in midwinter when basin yield totals exceeded total precipitation inputs. There was more net recharge in warm years than in cool years for most months in all regions, particularly for the SLP which had on average 19 mm more net recharge in warm years. Across the entire state, as well as each region individually, the 28 distribution of basin net recharge values between year types is statistically significant at the 95% confidence level for both tests (Table 2.4). Figure 2.6. Regional monthly net recharge and cumulative net recharge in different year types. Cumulative median net recharge (increase in storage) for each region and year type. By region: a) is the southern LP, b) is the northern LP and c) is the UP. In the LP warm years have more recharge earlier in the spring, and more total recharge at the end of the season. In the UP, spring recharge is significantly higher in cool years. The first day of water year for each month is listed below it in parenthesis on the x-axis. Unlike differences in snowpack, during the deeper winter months the UP shows the least difference between year types, while the net recharge in the Lower Peninsula can be centimeters more in warm years. Cool year net recharge for April within both the NLP and UP increases substantially in response to more melt from the persistent snowpack. This is reflected in the region’s very high S/P ratio in cool years (Table 2.3). Peak monthly net recharge amounts give insight to recharge mechanics in the different regions, with the peak monthly net recharge in the 29 UP making up the majority of the season’s total net recharge (79% in cool years and 57% in warm years), while the peak monthly contribution to net recharge decreases southward. In the NLP, the peak monthly net recharge contributed to 64% of the seasonal total net recharge in cool years and 48% in warm years; for the SLP it was 31% in cool years and 32% in warm years. Regional recharge amounts may also be influenced by the amount of liquid precipitation, with higher total net recharge in the SLP correlating with higher rain amounts. Net recharge in the Southern Northern Lower Upper Region Lower Peninsula Peninsula Peninsula Year Type Cool Warm Cool Warm Cool Warm Total Net Recharge 139 158 187 188 163 79 (mm) Peak Monthly 43 51 120 98 128 45 Recharge (mm) Total Rain (mm) 186 233 150 197 53 69 Mean S/P 0.42 0.30 0.58 0.46 0.83 0.68 Table 2.3. Regional recharge statistics. Statistics of estimated net groundwater recharge and precipitation for each region and year type. Cool years have more recharge in the LP, while the opposite is true for the UP. Warm years also have more rain in all regions and a lower ratio of snow to total precipitation (S/P). NLP, similar to streamflow differences, has less clear year type differences with total recharge similar in both warm and cool years despite differences in total rain and S/P, indicating the surficial geology’s role affecting recharge and not just melt/rain amounts. The deep surficial deposits with very high sand content leading to relatively high soil hydraulic conductivity values (2 to 10 mm/day for much of the NLP) could contribute to the lack of difference seen in streamflow and net recharge in different year types (Figure 2.1). Net recharge is higher in the Lower Peninsula than in the UP for much of the year, in all year types. This is likely related to more liquid winter precipitation in the southern regions. 30 Region Variable All SLP NLP UP Total Melt <0.001 <0.001 <0.001 0.2 Total Basin Yield <0.001 0.04 0.3 0.008 Peak Basin Yield <0.001 <0.001 0.6 <0.001 Day of Water Year BY50 Day of Water <0.001 <0.001 0.3 <0.001 Year Total Net Recharge 0.005 <0.001 0.02 0.02 Max SWE <0.001 <0.001 <0.001 <0.001 Max SWE Day of <0.001 <0.001 0.9 <0.001 Water Year Total Snow <0.001 <0.001 <0.001 <0.001 Peak Melt <0.001 <0.001 0.8 0.04 Peak Melt Day of 0.03 0.1 0.3 <0.001 Water Year SM50 Day of Water <0.001 <0.001 0.3 <0.001 Year Bare Ground Days <0.001 <0.001 <0.001 <0.001 Season Length <0.001 <0.001 <0.001 0.08 Table 2.4. Regional statistical analysis results. P-values from the Mann-Whitney ranked sum tests of difference on annual basin statistics grouped by year type. Italicized cells are p-values that are above the 95% confidence threshold for statistical significance. Kolmogorov-Smirnov analysis results always had p-values that were the same or lower than the Mann-Whitney analysis, so they were not included. 4. Discussion Our analyses support all four study hypotheses: 1) snow melt occurred earlier, and in lower quantities, in all regions during warm years; 2) complete melt events occurred more frequently in warmer years, leading to more days without snow on the ground; 3) basin yield (and thus stream discharge) peaked earlier and lower during warm year winters and finally; 4) there is more net groundwater recharge during warm years in the southern LP, but not in the UP. We also found that differences between warm and cool years were highest in the North and decreasing to the South. Average snowpack SWE peaks lower in warm years, with earlier melting of the snowpack (Figure 2.3). This difference between SWE peaks across year types is more notable in 31 the northern regions, likely due to the greater snowfall amounts these areas receive. Warmer years also produce less snowmelt, and thus a smaller snowmelt component of watershed hydrologic fluxes (Figure 2.4). How and when this melting occurs is also noticeably different between year types. Warm and cool years had similar numbers of melt events in all regions, but more of those melts are progressing to completion in warm years. The average melt event length and amount are less in warm years than cool years, especially in the UP. So while there may be similar numbers of melt events in both year types, the snowpack is less persistent throughout the season in warm years and melts more quickly. Indeed, for the UP there are 40 more bare ground days in warm years, indicating substantial portions of the season without snow cover, corresponding with increased melt events and complete melts. This differs from recent findings of Musselman et al. (2017) who indicated that warmer climate produces slower melt rates because more melt is occurring earlier in the year when days are shorter and the solar declination is lower. However, that study focused on thicker snowpacks in the western United States located at higher elevations which persist later into the spring, and as a result may be more sensitive to such climatic changes. Another facet of the changing snowpack dynamics is the response of Lake Effect snow amounts to climate warming. Burnett et al. (2003) found that warming temperatures in response to global climate change may be driving higher Lake Effect snow totals as the warmer temperatures lead to less ice cover and higher evaporation rates. Even with possibly increased amounts of evaporation from the Great Lakes under warmer temperatures, those same warmer temperatures are likely to lead to more of that Lake Effect precipitation occurring as rain rather than snow during the winter months (Champagne et al., 2020; Hayhoe et al., 2010). While snow amounts may be increasing through the decades, this study’s short temporal window is unlikely 32 to capture such climatic influences. Furthermore, the complex interactions of atmospheric oscillation indices, rising global temperatures and local variations in wind, humidity and other factors contributing to the Lake Effect phenomena is beyond the scope of this study. Streamflow, especially in the northern regions, responds dramatically to changing snowmelt in warm years. Streamflow patterns are responsive to earlier melting in warm years. Basin yields were lower and peaked much earlier in warm years, correlating with the earlier peak melt events in these years. These results also agree with what many studies of changes to snowmelt hydrology under warming climate scenarios have found (Boyer et al., 2010; Campbell et al., 2011; Hodgkins et al., 2003; Hodgkins and Dudley, 2006b; Johnston and Shmagin, 2008; Novotny and Stefan, 2007; Stewart et al., 2004). The influence of landcover on basin yield should also be considered in further studies to help elucidate regional differences in surface flows. For instance, the more urbanized southern regions may show less difference between warm and cool year flows because no matter when melting occurs, much of that melt will become overland flow on the more extensive impermeable surfaces. However, the role of land cover in these results is likely limited, as urban areas are generally limited (except in the SLP) and tend to be located lower within watersheds with most differences attributed to latitudinal climate gradients including significantly lower annual snowfall in the south, which leads to more intermittent snow cover. This also agrees with the results from Huntington et al. (2004), who found that the more northern New England study sites are likely to experience the most significant downward trends in snow to precipitation ratio (S/P) and earlier peak flows. The downstream effects of these streamflow changes are not entirely understood, and lake level differences between years should be further examined in future research since streams contribute close to half of the water stored in the lakes, and lake 33 levels have been projected to decline with warmer temperatures (Angel and Kunkel, 2010; Mortsch et al., 2000). Net recharge to shallow groundwater peaked earlier in warm years and exceeded the cool year net recharge in the southern regions. In warm years, the highest amount of net recharge occurs in March, a month earlier than in cool years. Indeed, the effect may be greater than calculated here, because some plants are already starting to exit their dormant phase in April, particularly in warm years. While we have assumed that winter ET rates are negligible for Equation 1, significant ET during April of warm years would reduce net recharge and make the shift in recharge timing even more pronounced. Kirchner and Allen (2020) determined from end- member “splitting” analysis of experimental data collected at the Hubbard Brooks Experimental Forest in New England estimated that little to no ET originated from snow season (Dec-Mar) precipitation (15 ± 15%). Using NLDAS-2 forcing data, we calculated average potential evapotranspiration (PET) for Michigan from November-April for the study years, finding the winter PET for the majority of the state totals around only 10 mm (Xia, Y., 2012; Figure A.2.2). Regardless of ET effects, by the end of April warm years still show higher overall net recharge than cool years in the southern LP, tentatively confirming the hypothesis of increased recharge in warm years. However, the hydrologic pathways for recharge to shallow groundwater are complex, and thus process-based modeling will be needed to fully evaluate these processes. The cause of this increased recharge is currently unclear and is further complicated due to the influence of frozen soil on infiltration rates. Several studies have shown that reduced snow depth and days with snow on the ground decrease thermal insulation of soil moisture, leading to more frozen ground and reduced infiltration into the soil (Isard and Schaetzl, 1998; Iwata et al., 2011). 34 Despite the increased likelihood of frozen soil in warmer years, there appears to be more net recharge in such years in the south. This could be due to preferential flow paths through and around frozen soils developing like those found by Mohammed et al. (2019) in the Canadian prairies, or it could be due to climatological and geological differences between regions. Regional differences similar to those found in the melt and basin yield results are present in net recharge estimates as well, with the southernmost region showing the least difference between year types. A number of potential mechanisms may produce these regional differences. For example, the general hydrologic processes governing runoff generation may be affecting how much melt infiltrates the subsurface instead of contributing to surface flows. Shallower snow depths that intermittently melt may allow proportionally more percolation to the water table compared to large melts that cause soil saturation early in the melt leading to overland flow. This may be part of the reason that the NLP had little difference in peak streamflow and center of volume arrival time; soils there are fairly uniformly sandy and coarse-textured with high hydraulic conductivity, compared to those of the SLP and UP (Figure 2.1; Soil Survey Staff, 2020). Mean hydraulic conductivity for the NLP is 2200 mm/day, compared to an average of 1000 mm/day for both the SLP and UP soils. Increased rain-on-snow events may also increase recharge as the warmer rain could increase the soil moisture temperature. In addition, more melt occurs in the late winter and early spring months in warmer years long before plants are active. Regional differences in melt and basin yield amounts also likely affect recharge. Again, due to the complex interplay of physical processes, process-based integrated hydrologic modeling of recharge is necessary improve the understanding of such relationships. Regardless of the drivers of these regional differences, there appears to be a latitudinal transition zone across the state, with the regions north of the transition zone (the UP) showing more 35 pronounced differences in snow hydrology between year types than the southern regions. The exception is total annual melt, for which all regions showed less melt (and correspondingly less snowfall) during warm years. The lack of difference between streamflow and recharge during warm and cool years in the south is likely because these regions receive less snow than the northern parts of the state. When examining the distributions of year type datasets regionally, the NLP appears to be the transition zone as it is the only region to consistently fall below the 95% confidence interval using both statistical tests (Table 2.4). This transition zone may represent a latitudinal gradient south of which receives too little annual snow for large timing shifts and north of which warming has very distinct melt effects. Thus, changes to snowmelt hydrology resulting from a warming climate aren’t as impactful to the overall water budget in these transitional regions. The more northern region showing a larger difference between year types is similar to findings by Suriano et al. (2019) who found that from 1960-2009 snow depth amounts across the Great Lakes Basin declined by approximately 25%, with the most significant decreases in the northern areas of the basin. An alternative explanation is that the global climate hasn’t yet warmed sufficiently to impact the hydrology of the northern regions in cool winters as it has the south. To confirm this, study of longer-term trends is needed to establish the location and stationarity of this transition as the climate continues to warm. 5. Conclusions Using observations from GHCN and the USGS, alongside model reanalyses from SNODAS and PRISM, this study found that Michigan’s snowmelt hydrology shifts significantly between years defined as “warm” and “cool”, with the degree of shift varying by region. In all regions, warm years had less total snowmelt with earlier melt occurring in more complete events than cool years. As expected, precipitation also shifted towards more rain in warm years. These 36 changes to melt and precipitation caused earlier and lower peak streamflows in warm years, with less groundwater recharge. Differences between regions indicate a “transition zone” where southern regions show smaller differences between warm and cool years. Future studies will focus on expanding the spatial and temporal scope of this research. The most significant limitation in this study was the short timeframe available from the snow model. To look further back in time, the only source of empirical climate data is historical weather data, which has spatial limitations and only provides snowfall/snow depth data without the SWE component. Historical climate model outputs are available, but they either don’t include the SWE component or are too coarse in their spatial resolution. By using the metrics defined here, the multimetric analysis can be applied to years over a longer period to robustly classify warm versus cool years. Then using this multimetric analysis with stream gauge records for those expanded years, inferences about melt processes can be made without having actual melt data—which are generally poorly available prior to SNODAS. We plan to apply this process across larger scales to better understand how widespread these hydrologic changes may be. Ultimately the results of this and future work will form an observational foundation for targeted and improved simulations of these melt processes to quantify how winter and spring hydrology in the Great Lakes will likely change in the coming decades due to global change. Acknowledgments The text of this chapter is a reprint of the material as it appears in Journal of Hydrology in a manuscript first published September 10, 2020 (doi:10.1016/j.jhydrol.2020.125517) coauthored by David Hyndman and Anthony Kendall. This research was supported by USDA- NIFA grant 2015-68007-2313, NOAA grant NA12OAR4320071, and NASA grant NNX11AC72G along with Michigan State University’s Environmental Science and Public 37 Policy Program and the Department of Earth and Environmental Sciences. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of USDA NIFA, NASA, or NOAA. 38 APPENDIX 39 Figure A.2.1. Normalized seasonal temperature across the domain. Time-series of the mean annual metric score compared to the average daily temperature and minimum daily temperature across the entire state for three different periods: the entire snow season (Oct-May), just the winter (Dec-Feb) and just spring (Mar-May). Temperature values are normalized to the 13-year mean. 40 Figure A.2.2. Winter PET for study region. Average November-April potential evapotranspiration (PET) for years 2003-2017. Other than near-shore environments the majority of the region’s PET during the winter months only totals between 10-20 mm. Data for this figure comes from NLDAS-2 forcing data (Xia, Y. et al., 2012). 41 CHAPTER 3: SNOWPACKS DECREASE AND STREAMFLOWS SHIFT ACROSS THE EASTERN US AS WINTERS WARM Abstract Climate change is increasing winter temperatures across the planet, altering snowmelt hydrology. This study addresses a gap in snow research in non-alpine areas by examining changes to snow and winter and spring streamflow across most of the eastern US using daily observations from weather stations and stream gages from water years 1960-2019. These daily data were aggregated across drainage basins and classified winters with similar temperatures; differences between winters and both seasonal and annual trends were statistically quantified. Winters were classified as “warm” or “cool” in each drainage basin relative to the 60-year mean; analysis of the data indicates that warm winters occur more frequently in recent decades from an average of 0.39 to 3.96 warm winters/decade from the 1960’s to the 2010’s respectively. Those classifications were then used to examine changes in snowpack over the same period, which shows that warmer winters have on average 50.1 cm less annual snowfall, a reduced maximum snowpack depth by 14.4 cm, and 34 more bare ground days. These changes correlate with shifts to higher winter streamflows as well as peak basin yields that are 0.02 cm lower and occur three days earlier in warm winters. In addition to altered soil moisture and stream ecosystem dynamics, these snow and streamflow changes may have negative infrastructure and economic implications including impacts to winter tourism and agriculture. 42 1. Introduction Warming winter temperatures are occurring around the planet driven by climate change in response to anthropogenic activities. These warmer winters will affect snowmelt regimes, with significant implications for winter and spring hydrology in areas where snow dominates winter precipitation (Adam et al., 2009; Hyndman, 2014). Numerous studies have already documented these climate-driven snow regime changes, but much of that research has focused on arctic or alpine regions of the planet where snow accounts for much of the annual water budget. However, the expansive mid-latitude non-alpine regions of the globe are also seeing changes to perennial snowpacks and have been relatively understudied. These areas, often with abundant shallow groundwater, are responsive to changes in snowmelt dynamics since the spring snowmelt is a large driver of peaks in annual streamflow. Thus, it is critical to more fully examine these understudied areas to better understand associated hydrologic changes, which are important for stream ecosystem health, water resources, flood mitigation, and hydrologic infrastructure design. Changes to alpine snowmelt and associated streamflow due to warming temperatures have been documented in numerous studies of sites across the world. Most of these find similar trends in snow including declining snowpack thickness (Hamlet et al., 2005; Mote, 2003), earlier snowpack melting (Clow, 2010; Hamlet et al., 2005; Ishida et al., 2019; Jefferson et al., 2008; Stewart et al., 2004), and precipitation shifts to more rain and less snow (Ishida et al., 2019; Mote, 2003). These changes in the snow dynamics have led to earlier seasonal flows (Clow, 2010; Jefferson et al., 2008), earlier peak flows (Jefferson et al., 2008; Stewart et al., 2004), and earlier center of discharge volume (CDV) timing (Hidalgo et al., 2009). Climate projections generally agree that a warmer climate in the mountainous regions of the globe will exacerbate these changes through decreased snow cover and melt runoff (Javadinejad et al., 2020), earlier 43 melt (Brubaker and Rango, 1996) and less contribution of snowmelt to streamflow (Li et al., 2017). Many studies have examined changes to snow at broad spatial scales and have found similar results to those from alpine settings. Warmer winter temperatures across North America are leading to reduced snow cover (Dyer and Mote, 2006; McCabe and Wolock, 2010; Gan et al., 2013) and lower snow/total precipitation (S/P) ratios (Feng and Hu, 2007; Berghuijs et al., 2014). Such changes in the Northern Hemisphere are projected to continue under a warming climate with most models showing decreased snow cover and melt (Adam et al., 2009; Barnett et al., 2005; Peacock, 2012), continuing decreases in S/P (Mankin and Diffenbaugh, 2015), and decreased late-season snow (Manking and Diffenbaugh, 2015) leading to earlier peak streamflows (Barnett et al., 2005; Clow, 2010). Much of the central and eastern United States receives substantial snowfall, with a historically persistent snowpack. Several studies have examined snow and streamflow dynamics in this region. In New England and the northeastern US there have been more bare ground winter days (Burakowski et al., 2008), decreased snowfall amounts (Burakowski et al, 2008), decreased snow depth (Hodgkins and Dudley, 2006a; Hayhoe et al., 2007), decreased S/P (Huntington et al., 2004; Javed et al., 2019), earlier melting (Campbell et al., 2011) and declining snowpacks (Campbell et al., 2011). These factors have been contributing to earlier and reduced peak flows (Campbell et al., 2011; Hayhoe et al., 2007) and earlier seasonal flows (Hogdkins and Dudley, 2006b). Similar observations of changes to snowmelt and streamflow have been made for the Great Lakes region (Ford et al., 2020; Hayhoe et al., 2010; Hodgkins et al., 2007; Johnson and Stefan, 2006; Suriano and Leathers, 2017; Suriano et al., 2019). 44 The future of snow hydrology in the eastern US under a warmer climate has similar future projections as alpine settings. In some regions of the eastern US, climate change is expected to increase winter precipitation (Byun and Hamlet, 2018; Naz et al., 2016), but more of that precipitation will occur as rain rather than snow (Boyer et al., 2010; Byun and Hamlet, 2018; Chin et al., 2018; Hayhoe et al., 2010). The snow cover will be reduced (Demaria et al., 2016; Hayhoe et al., 2007; Mortsch et al., 2000) as will the snow-water equivalence (SWE) and melt (Pradhanang et al., 2011; Naz et al., 2016) and that melt will occur earlier in the season (Tu, 2009). These snowmelt changes will lead to reduced spring runoff from melt (Mortsch et al., 2000), earlier seasonal flows (Champagne et al., 2019; Tu, 2009), earlier CDV time (Boyer et al., 2010), earlier peak flows (Naz et al., 2016) and possible increases in extremely high and low flows in winter and spring (Naz et al., 2016). Here changes to snow and streamflow are examined across the northern half of the central and eastern Continental United States. Specifically, historical weather station and stream gauge observations from water years 1960-2019 are analyzed within drainage basins in 32 states ranging from the eastern Atlantic Ocean coast to parts of the Dakotas and Nebraska and from the northern US-Canadian border to the fringes of the American South. These sites were chosen to encompass the regions of the eastern US that receive enough annual snowfall to possibly have measurable effects on spring streamflow. Despite the numerous studies on climate change impacts on snowmelt hydrology, there is limited research published on such changes in non-alpine settings. Most of the studies of the non- alpine settings in the United States have either focused on subregions (e.g., New England, Great Lakes, Midwest) or have only focused on one component of the system (i.e., only snow or streamflow). An exception was the study by Hodgkins and Dudley (2006b) that examined both 45 snow and streamflow changes across most of the eastern United States, but the observational period for the study ended in 2002, which is missing the most recent decades during which there has been significant climate change relative to long-term norms. While there has been limited scientific research into the climate-driven effects on snowmelt hydrology in these non-alpine settings, these areas will likely experience significant environmental changes as many drainage basins in this area receive significant contributions to their annual water budget from snowmelt. This paper addresses a gap in the literature by quantifying changes to snowfall, snowpack depth and streamflow by classifying winters as warm or cool and examining the differences in snowmelt hydrology between these winter types. The hypotheses examined in this study are: 1) winters with warmer temperatures are becoming more common, causing 2) decreased snowfall and shallower snowpacks, 3) leading to more days with bare ground and changes to streamflow in the winter and spring, such as 4) higher winter flows and 5) earlier and reduced peak flows. These hypotheses are examined using 60 years of daily observational data of temperature, snow and streamflow to identify warmer winters and examine hydrologic changes in those seasons compared to cooler winters. The manuscript starts with a description of the relevant climate, hydrology, and geology of the study region. Then, the stream gauge and climate station data sources used are described along with the classification, aggregation, statistical and trend analyses conducted on that data. The decreased snowpack depths, earlier melting, increased bare ground and earlier streamflow trends found by these analyses are then presented, followed by a discussion of the implications of the results that are contextualized within the existing literature. Finally, the research findings are summarized, and the next steps of simulation modeling are described in the conclusions. 46 2. Methods 2.1 Study Region Figure 3.1. Study region maps. a) Mean annual snowfall amounts from GHCN stations averaged across Hydrologic Unit Code 4-digit (HUC-4) basins used in this study. Average snowfall amounts range from less than 35 cm per year to over 300 cm per year. b) shows the HUC-2 basins used for grouping in this study, as well as the 1,725 GHCN stations that provided temperature and snow data. The eastern and central (hereafter just “eastern”) United States has a diverse range of geology, land use, and climate conditions. This study region stretches from the sandy prairies in eastern Nebraska and Kansas to the mid-Atlantic coast, from the clay-rich fertile soils of the agricultural interior to the heavily glaciated lands of the Great Lakes (Figure 3.1). These boundaries were selected to capture areas with enough annual snowfall to contribute 47 substantially to the annual water budget, but otherwise encompass generally non-alpine geography. The study region includes parts of the Appalachian Mountains, the five Great Lakes, and the agricultural lands of the Mississippi drainage basin. The area is geologically diverse with thick, clay-rich surficial deposits in the southern Midwest that transition to sandier soils in the Great Lakes before becoming thin soils underlain by Precambrian igneous rocks of the Canadian Shield. The eastern boundary of the region is the Atlantic Ocean, capturing the Appalachians. While the Appalachians are topographically distinct from the surrounding region, they are relatively subdued compared to the younger, taller mountains of the western Cordillera and the orographic effect in this eastern portion of the study region is not as consequential. The land use varies as much as the geology across the region. The southern and western parts of the study region are mostly large-scale row-crop agricultural fields interspersed with tracts of deciduous forests. The southern Great Lakes region are similarly heavily agricultural, with the northern regions more dominated by coniferous forests and wetlands. Much of the urbanized land in the region is in the southern and eastern sections, potentially affecting stream response to melt changes as more impervious surfaces generate more overland runoff. The climate across the study area is generally humid, receiving consistent precipitation in most months of the year. Annual snowfall amounts increase northward from the southern portion of the study area that receives less than 35 cm/year to some northern portions receiving upwards of 300 cm/year (Figure 3.1). Across the region, mean annual precipitation is 82.6 cm, with 49.5 cm of that falling during October through May (Figure A.3.2). Annual snowfall amounts in the northern areas are also regional, with the highest snowfall amounts occurring in the Great Lakes and the Northeastern states. High snowfall amounts in the Great Lakes states are driven by the 48 Lake Effect phenomena in which evaporation from the lakes contributes to larger snow amounts on the downwind side of the lakes. High snowfall events occurring in the Northeastern US are typically caused by cold Arctic air interacting with the warmer Atlantic Ocean waters, generating the cyclonic storm systems regionally known as “Nor’easters.” Temperatures vary seasonally and regionally, but generally remain near freezing for much of the snow months (Supplemental Figure 3.3). Mean winter (DJF) daily average temperatures across all sites and years is -3.5 °C. Mean spring (MAM) daily average temperature is 9.1 °C, but mean spring daily minimum temperature is only 2.7 °C, with the northern stations at or below freezing in the spring. Across all sites and years average days above freezing in winter months is 31 days, while average days above freezing for the snow season (Oct-May) is 162 days. Streamflow in most the study area is relatively consistent throughout the year. Major stream systems in the region are perennial, with many of the rivers sustaining flows in months with low precipitation amounts due to shallow groundwater discharge in the streambed. Annual peak flows typically occur in the spring from large inputs due to snow melt. Flows then recede with periodic pulses in response to precipitation events in the summer and early fall before stabilizing at baseflow levels for the rest of the year. Streamflow is typically responsive to precipitation events, depending on local soils, land use and topography. 2.2 Data Sources All data used in this study were derived from daily observational data provided by either the US Geological Survey (USGS) or the Global Historical Climatology Network (GHCN). Daily station values of temperature, snowfall and snow depth were downloaded from the GHCN station records using the rnoaa software package (Chamberlain et al., 2020), while daily values of average temperature were calculated from daily maximum and minimum temperature station 49 data. Streamflow data for this study are compiled from daily observations of stream discharge across the USGS gage network using the software package dataRetrieval (DeCicco et al., 2020) (USGS, 2018). These daily observations were collected from October 1, 1959 – May 31, 2019 (water years 1960-2019, here defined as October 1 of one year through September 30 of the following). This study focuses on the October-May period of each year because this period captures the significant snow accumulations in the region in all but the most extreme winters. Data were quality checked for both erroneously high values and continuity. Temperature data were removed if daily minimum temperatures were higher than the daily maximum temperature. Each station’s annual snow season (here defined as October-May) data were analyzed, and years missing more than 5% of the days during that period were removed. Daily snow data with amounts higher than 250 cm were considered too high to be accurate in this generally low-relief study region and removed based on the distribution of snow data values. This left data from 1,369 GHCN temperature stations, 1,725 GHCN snow and precipitation stations, and 1,751 USGS streamflow gauges for use in this study (see Figure A.3.1). 2.3 Analysis Methods Following data synthesis and QA/QC, the data were reduced via spatial and temporal aggregation prior to conducting further analyses and statistical tests. 2.3.1 Spatial and Temporal Aggregation for Analysis Three distinct spatial scales were used for the analysis: 1) all stations over the study region, 2) summaries within Hydrologic Unit Code 2-digit (HUC-2) polygons (Seaber et al., 1987), and 3) summaries within Hydrologic Unit Code 4-digit (HUC-4) polygons (Seaber et al., 1987). The station/gauge point data were imported into the “R” (R Core Team, 2020) 50 programming software, then combined with rasterized shapefile polygons of the HUC-2 and HUC-4 basins to extract basin averages. For each of the HUC-2 and HUC-4 scales, station values were averaged annually within polygons prior to other analyses. Note that USGS gauges have distinct basins from those delineated by the HUC system, however those gauges were treated for this analysis as lying entirely within the HUC-2 or HUC-4 basin. In general, this is a good assumption as these large-scale basins are much larger than individual stream gauge basins. Four temporal aggregations were used for these analyses: 1) all years individually, fit to linear regressions (described below), 2) “warm” and “cool” winters as classified via multimetric score, 3) decadally-binned values, for all years as well as the warm and cool subsets, and 4) monthly-binned values were used in the analysis of the precipitation data. These temporal scales will be identified both in the methods described below as well as in results section. Figure A.3.4 provides a visualization of the methods in a flow chart. 2.3.2 Multimetric “Warm” and “Cool” Snow Season Classification Snow seasons were classified as warm or cool for each basin using a multimetric analysis that combines six temperature related metrics: average and minimum temperature averaged for the entire snow season (Oct-May), minimum and average temperature averaged for just the winter months (Dec-Feb), and the spring months (Mar-May). These metrics were calculated for every HUC-4 watershed in each year, and then normalized to the 60-year median for each watershed; more positive scores thus indicate temperatures that were warmer than the long-term median and negative scores cooler. This normalization calculated the difference between the annual value and the 60-year median, then divided that difference by the 60-year standard deviation. Then, the arithmetic mean was taken across each of the six metrics for each year and basin giving each HUC-4 an annual score. For visualization, the arithmetic average score was 51 calculated for each year across the study domain as the arithmetic mean of the score across all GHCN stations in the study domain. Trends within the scores were calculated via linear regression, both at the basin scale and across all stations. For these regressions, the multimetric score was regressed against year, allowing for an intercept term. Annual scores greater than or equal to 0.5 were classified as “warm” and less than or equal to -0.5 as “cool”. This threshold was chosen because it indicates snow seasons that are warmer or cooler, while not being too heavily influenced by outlier precipitation events. Since the metric score was calculated for each basin, the classifications are basin-specific (so a specific year may be classified as warm or cool for some basins but not all). This classification method was reasonably accurate at identifying warmer than normal winters for data from a much smaller region and shorter temporal scale (Ford et al. 2020). Note that throughout this manuscript, references to warm and cool “winters” refer to the entire Oct-May snow season unless otherwise specified. Similarly, the phrase “winter type” refers to whether a snow season has been classified as either cool or warm according to this scheme. 2.3.3 Regional-Scale Snowpack, Streamflow and Precipitation Analyses Using the classifications from the multimetric analysis, snowfall and snow depth data from GHCN stations were analyzed within winter types. Daily snow data from GHCN stations were grouped by winter type, then statistics were calculated for annual snowfall amounts, peak snow depth, and bare ground days, which are defined as those with snow depth less than 2.54 cm. Winter type stats were then spatially aggregated across HUC-4 basins for visualization and mapping. Average day-of-year snow depth was also extracted across winter types and the larger HUC-2 basins (Figure 3.1) to examine regional differences. 52 Once snow differences between winter types had been analyzed, streamflow was examined through the warm and cool winter lens. Daily discharge values were first converted to daily basin yield by dividing the discharge by the basin drainage area so that flows from basins of different sizes could be more readily compared. Annual statistics were then calculated for each gauge site and compared across winter types. Streamflow data computed include peak basin yield timing and amount, along with the timing of the CDV of flow. To examine regional streamflow patterns in a similar fashion to snow depth, median basin yield across winter types within each HUC-2 was calculated from ten-day moving averages of daily basin yield for each gauge (to smooth out daily fluctuations in flow). Finally, the percentage of spring flow relative to the average flow of the rest of the year was computed for each HUC-2 basin through time to examine long-term trends in spring flow relative to non-spring flow. For each basin, the average basin yield in the spring months (MAM) was divided by the average basin yield in all other months (including non-snow months not included in other analyses in this study). For further context, these small values were then multiplied by 100 to convert the units to “%/year” and then multiplied by 10 for units in “%/decade”. GHCN precipitation data was also analyzed alongside the streamflow analyses results. Daily precipitation totals were first aggregated across HUC-2 basins and winter types. The aggregated daily data was then summed to monthly totals. 2.3.4 Distribution Analysis and Significance Tests The distributions of snow and streamflow data between winter types was examined using two statistical tests: the Mann-Whitney (Wilcoxon) ranked sum test (Hollander and Wolfe, 1973) and the two-sample Kolmogorov-Smirnov (KS) test (Conover, 1971). First, HUC-4 seasonal averages were calculated for snow and streamflow variables such as peak snowpack depth and 53 time, season length, peak spring flow amount and time, and total seasonal flow. These basin averages were then grouped by winter type and the dataset was checked for a normal distribution using the Shapiro-Wilk Normality Test (Royston, 1995). Data that were deemed a normal distribution in the results of the Shapiro-Wilk test were confirmed to be normal by plotting values in a histogram. Once the data were determined to be normal or non-normal the ranked sum test and KS test were applied to examine statistical differences between the two year-type datasets. Variables with normal data were also analyzed using the two-sided Student’s T-test. To determine if the datasets are different the p-values from all tests were evaluated on a 95% confidence threshold. 3. Results 3.1 Winter Type Classification and Trends Several spatial and temporal trends emerged from the multimetric winter type analysis. The number of warm years per decade is increasing (Figure 3.2a); most basins only have one or two warm winters in the 1960s, while multiple basins experienced eight warm winters in the 2010s. Not only does the number of warm winters increase across basins, but so does the average annual metric score within each decade (Figure 3.2b). Most basins in the 1960s and 1970s have average winter scores that are cooler than the longer-term mean, but by the 1990s most basins have warmer-than average winter scores. Certain regions show more warming than others, with the eastern-most basins in the study area having particularly warm winters in recent decades. 54 Figure 3.2. Regional warm winters and metric scores per decade. Maps of HUC-4 averaged a) warm winters per decade as defined by the multi-metric analysis, b) decade-averaged metric score for each basin. Noted in each map is the cross-basin average count or score for each decade. Note both the increasing trends in the number of warm winters and the mean metric scores per decade across all basins. 55 This positive temporal trend of warmer winters is apparent within the average multimetric score across the study domain (Figure 3.3). This is true when examining the decadal average increases as well as the linear trend across all years. It is worth noting that average annual scores appear to deviate further from the 60-year median through time in both winter types, with both the warmest winter on average (2012) and the coolest winter on average (2014) in the last decade. This is represented by the increased standard deviation away from the decadal median displayed as the grey-shaded area in Figure 3.3. Figure 3.3. Metric score time series. Mean winter metric score across the study domain. Each point represents the average of the multi-metric score across all basins for that year. The red and blue lines represent the threshold for warm and cool winter classifications, respectively (though elsewhere winters are classified as warm or cool at the basin level, rather than across the entire domain). The green line represents the decadal average. The black line is the linear trend of the data for all years. The gray shaded squares represent the standard deviation of scores for that decade. The average increases through time, as does the spread of the data with more warmest and coolest scores occurring in the more recent decades. 56 3.2 Snow Depth and Precipitation Analysis Warmer winters generally have less snow, which covers the ground for fewer days and melts earlier than in years with cool winters. Annual snowfall totals in warm-winter years were 50.1 cm lower than in cool-winter years across all basins (Table A.3.1). The reduced snowfall amount likely contributes to a thinner snowpack throughout the snow season, with warm winters having a 14.4 cm lower peak snow depth than in cool winters, with that peak occurring 12 days earlier. This thinner snowpack is also less persistent, with warm snow seasons having 202 bare ground days compared to 168 days in cool snow seasons across all basins. The season length, here defined as the amount of time between the first day of snow >= 2.54 cm and the last day of snow >= 2.54 cm, is also shorter by 24 days in warm winters. When the Shapiro-Wilkes test was applied to the data, only the peak snowfall day of water year data was confirmed as having a normal distribution. Using the 95% confidence threshold for the statistical tests of difference, all snow variables had significantly different distributions between winter types (Table A.3.2). 57 Figure 3.4. Regional bare ground days maps. The average amount of bare ground days in a year (Oct-May) for each HUC-4 basin in a) the first decade of the study period (1960’s) and b) the last (2010’s). There is a northward progression of bare ground days as more northern regions have more bare ground days per winter through time. Panel c) shows the slope of the bare ground days per season through time for each HUC-4, with most showing an increasing amount of bare ground days. Out of 84 basins, only 2 exhibited a decreasing slope in bare ground days per season. There is a northward migration of the number of bare ground days through time, with the northern regions having more bare ground in recent decades (Figure 3.4). Basins had 27 more bare ground days in the 2010s than in the 1960s. This is represented by the northward advancement of basins with 160 or fewer bare ground days (which corresponds to 80 snow days during the 240-day snow season period) in Figure 3.4. In some areas this boundary migrated north by several hundred kilometers by the end of the five-decade study period. 58 Bare ground days per snow season are increasing in nearly all basins, visualized by the slope of bare ground days per season (Figure 3.4c). Changes are occurring most rapidly across a wide central swath of the study region; slopes are largely uniform longitudinally. There is an increase in number of bare ground days through time in regions in the interior of the continent (e.g., the plains states, the Midwest and much of the Great Lakes area), while some of the southernmost and northernmost basins in the region show little increase to even slight decreases. Notably, patterns in the slope of bare ground days per season are largely uncorrelated (r2 = 0.31) with a similarly calculated (Figure A.3.5) slope in multimetric score per season. The rate of change in bare ground days is highest where average winter temperatures are closest to 0 in recent decades. Typical snowpack depths show dramatic differences in both thickness and duration between warm and cool winters across the study domain (Figure 3.5). As expected, these regional winter type averages of daily depth show earlier melting of the snowpack in the south, for both warm and cool years. Only 7 of the 9 HUC-2 regions had reliable annual snowpacks; the southernmost (Arkansas-White-Red and South Atlantic) experience episodic snow events in both warm and cool winters. The seasonal snowpack melts substantially earlier in warm winters for those 7 HUC-2 regions with reliable annual snowpacks. Differences between typical warm and cool winters range from 16 days in the New England HUC-2 to 41 days in the Missouri HUC-2. Similarly to melt dates, median daily snow depth in warm years across all HUC-2 regions is lower, with greater differences between winter types for the more northern regions. Average maximum basin snow depth decreased by 50% or more in two HUC-2 basins, with eight out of nine basins decreasing by 30% or more. Peak snowpack depths also occur much earlier for warm years, ranging from 1 to 29 days. 59 Across basins, precipitation in different winter types varied little. Many basins showed slightly more precipitation in the spring months (MAM) in warm winters (Figure A.3.3). The pattern of monthly precipitation for the rest of the year was similar for most year types with the lowest amounts of monthly precipitation in the winter months (DJF) and the highest early in the summer. Annual totals of precipitation were slightly higher in warm winters in all but one basin, but the difference between annual precipitation totals of different winter types was 5 cm or less in all the basins. Figure 3.5. Regional winter type snow depths. Mean daily snow depth averaged across winter types and HUC-2 basins. The panels are roughly arranged by latitude of the HUC-2. The basins are: a) Souris-Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic. Note the different y- axis scales for each row. A small inset of the HUC-2’s with that panel’s basin filled in is included for spatial referencing of the data. Warm year depths are noticeably lower for most seasons in all but the southern-most basins. For reference, Day 1 of water year is October 1st and day 250 is June 7th. 60 3.3 Streamflow Analysis Streamflow occurred earlier in warm winters with decreased peak flow amounts and increased annual flow totals compared to cool winters. Across the entire study region, during warm winters the median peak basin yield occurred 3 days earlier and was over 14% lower (0.02 cm/day) than in cool years (Figure 3.6; Table A.3.1). The arrival of half of the annual runoff volume, measured by the CDV statistic, occurred earlier in warm winters for all HUC-2’s except the Arkansas-White-Red with a median CDV across basins occurring 14 days earlier (Figure 3.7). Median total basin yield amounts increased by 13.7% (2.9 cm/day) in warm winters compared to cool when taken from station annual totals. All streamflow variables were non- normal when assessed using the Shapiro-Wilkes test. Using the 95% confidence threshold for the statistical tests of difference, only one variable was not significantly different between winter types: Seasonal Basin Yield Amount (Table A.3.2). These statistical tests were applied to station data aggregated across basins rather than just the station annual values, which is likely the cause of the discrepancy between the statistical test results and the 13.7% increase in warm winter total basin yield. All other variables tested had significantly different distributions of warm and cool winter data according to the results of each of the tests applied. 61 Figure 3.6. Regional winter type basin yield. Each panel shows the 10-day moving average of daily basin yield averaged across winter types and HUC-2 basins. The panels are arranged roughly by latitude of the HUC-2 basin. The basin names are: a) Souris-Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic. A small inset of the HUC-2’s with that panel’s basin filled in is included for spatial referencing of the data. Note the different y-axis scales for each row. Most basins show higher flows earlier in the season in warm years but higher peark flows in cool years. Regional patterns emerge from examination of flows within basins. The more northern and eastern basins have higher flows regardless of winter type compared to the southern and western basins (Figure 3.6). The New England, Ohio and Mid-Atlantic basins have peak flows that are 0.052-0.096 cm/day higher and occur later in cool winters, in some cases by several weeks. The Souris-Red-Rainey and Great Lakes basins also have higher and later peak flows in cool winters, but not by as much (0.031 and 0.021 cm/day higher, respectively). The peaks in streamflow in these northern basins correspond to the timing of significant snowpack melt, especially in cool years (Figure 3.5). The remaining basins show similar patterns but the difference between winter types is minimal, especially in the southernmost basins with very low 62 basin yield amounts like the Missouri and Arkansas-White-Red (0.001-0.009 cm/day higher in cool years). Figure 3.7. Regional winter type cumulative basin yield. Cumulative basin yield averaged across winter types and HUC-2 basins. Panels are approximately arranged by latitude of the HUC-2 for the: a) Souris-Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic Basins. A small inset of the HUC-2’s with that panel’s basin filled in is included for spatial referencing of the data. Note the different y-axis scales for each row. Most basins show higher winter flows in warm years, with cool year flows not reaching the cumulative warm year flows until April when the spring melt occurs. The timing differences between winter types are even more apparent when viewed cumulatively; warm year flows shift earlier with more flow occurring in January and February (Figure 3.7). In six of nine basins cumulative flow in warm winters is higher than in cool starting 63 late in the winter months, and then stays higher the rest of the year (Figure A.3.7). Basin yield totals in these basins by the end of the water year is higher in warm winters, especially in the eastern basins. Median annual basin yield in warm winters is 2.6 cm higher in the Mid-Atlantic basin and 2.92 cm higher in the South Atlantic basin. The basins that do not follow these winter- type patterns like the Souris-Red-Rainey and Missouri generally have low basin yields regardless of year type. Figure 3.8. Regional spring vs annual discharge. Each panel is the scatterplot of the mean annual basin yield in the spring months (MAM) divided by the mean basin yield for all other months and then converted to a percent with a linear trendline fitted to the data. The basin names are: a) Souris-Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic. All basins except the two southernmost basins (g. and i.) have decreasing trends through time. 64 To better assess this timing shift in streamflow, the spring flow is compared to the rest of the year’s flow for each of the HUC-2 basins (Figure 3.8). All but the two southernmost basins, the Arkansas-White-Red and the South Atlantic, show decreasing linear trends through time with slopes ranging from -0.43 to -3.6 %/decade (Figure 3.8), all trends were significant with p-values < 0.001. In contrast the median flow in winter months compared to flow in the remaining months shows an increasing trend in all but the same two southern basins, the Arkansas-White-Red and the South Atlantic (Figure A.3.6). Slopes of winter flow through time for the remaining basins ranged from 0.08 to 0.80 %/decade with p-values for the trends in all basins again <0.001. 4. Discussion The results of the analyses confirm the study’s hypotheses. On a decadal scale, the occurrence and magnitude of warming winters has increased for most basins, especially recently compared to the 1960s and 1970s. Warmer winters typically have less snowfall and more frequent melting (shown also in just the Great Lakes HUC-2 over a shorter period in Ford et al., 2020), leading to decreased depths. The snow cover in these winters is less persistent throughout the season, covering the ground for fewer days and melting earlier. These snow changes correlate with differences in surface flows during warmer winters including shifts to earlier seasonal flows leading to higher winter (DJF) flows, lower spring (MAM) flows, and earlier center of discharge volume (CDV). Over the 60-year study period, the changes to snow because of warmer winters are clearly evident, and if warmer winters continue to become more common, shifts away from cool winter snow patterns will likely increase in most basins. As expected, the southern regions of the study see more bare ground days regardless of winter type, with most October-May periods lacking snow cover. These differences between winter types in the central and northern regions 65 is becoming larger as there is a northward shift of bands of similar bare ground days through time (Figure 3.4). The largest increases in bare ground days through time occur in a centrally located band of basins in the study region stretching from eastern Nebraska through the southern Great Lakes and extending to Mid-Atlantic states. This agrees with other studies of snow cover in the eastern US and at the continental scale (Burakowski et al., 2008; Dyer and Mote, 2006; Gan et al., 2013; McCabe and Wolock, 2010; Suriano et al. 2019). Suriano et al. (2009) found mean snow depth across the Great Lakes Basin decreased by 25% from 1960-2009, which is approximately twice the change found in this study for the same period (12%). However, the data in that study came from a gridded snow depth product instead of the raw station data used in this study. It should be noted that while this study assumes air temperature is the main driver of the snowmelt changes seen other factors are also influential in the role of snowmelt generation such as relative humidity and shortwave radiation (Ishida et al., 2019). This increase in bare ground days in warmer winters correlates with earlier melting of the snowpack. As with bare ground days, the southern basins show only slight differences between winter types, likely due to the low snowfall amounts regardless of winter type. The northern basins however show distinct differences in snowpack between winter types, with the warm winter snowpacks arriving later, melting earlier and having peak depths that are barely half the peak depth in cool winters. Regional studies of snow in the eastern US have found similar changes to snow including shallower snowpacks that melt earlier in the spring (Campbell et al., 2011; Ford et al., 2020; Hayhoe et al., 2007; Hodgkins and Dudley, 2006a Hogdkins et al., 2007; Johnson and Stefan, 2006). Hayhoe et al. (2007) showed the number of days with snow cover in winter months (DJF) in the northeastern US from 1960-1990 changed (decreased) at a rate of - 0.04 to -0.07 days/month/decade. Applying the same methods to this study’s data found a higher 66 rate of decrease of -0.13 days/month/decade for the same period. While the relationship between melt and surface flows are complex and often uncertain, changes to snow correlate with streamflow changes indicating a causal relationship. This shift to earlier melting in warmer winters is likely directly influencing the higher winter flows seen in these years. Figure 3.7 shows the cumulative basin yield averaged across HUC-2 basins and winter types. Most basins see warm winter flows shifting earlier, breaking away from the cool winter flows in mid-winter and staying higher until late in the season when cool winters have their peak melts in March and April. This shift can be seen in the 10-day moving average of basin yield taken across basins and winter types, where warm winter flows are higher for the first several months of the snow season, but the spring peak flow in cool winters occurs in April and rapidly rises above warm winter flows during that time (Figure 3.6). These findings are consistent with other studies of streamflow in the eastern US (Campbell et al., 2011; Ford et al., 2020; Hayhoe et al., 2007; Hodgkins and Dudley, 2006b). Hodgkins and Dudley (2006b) examined changes to streamflow timing in the northeastern US from 1913-2002, finding the CDV time from January-June advanced by 1-15 days from 1953-2002. Examining this study’s data for the same area from 1960-2002 revealed an advance of the CDV time from January-June by 7.5 days. These changes are also reflected in the amount of streamflow in the spring relative to the rest of the year, with most basins showing a decreasing trend through time (Figure 3.8). The two southernmost basins are the only ones to not show this decreasing trend because they receive very little annual snowfall relative to the more northern basins, thus this seems to indicate the role of melt in these changes to spring flows. It is unclear how these timing shifts will correlate with increased winter flooding and decreased spring flooding in warm winters. While peak flow is occurring earlier in warmer winters, those peaks are also decreasing 67 on average. The effect of melt changes on shallow groundwater will also partly determine winter flood prominence and needs to be more thoroughly investigated to better assess the risks. These analyses examined a facet of the basin-scale water balances driven by snowmelt. There are several critical portions of the water balance that were not included in this study but are influential to the observed changes. The role of land use and land cover change through time was not examined here, but 60 years is enough time for the landscape changes due to anthropogenic factors to play an important role. The soil freeze-thaw dynamics may also play a role, but this is difficult to fully quantify since the reduced snow cover reduces insulation of the soil from fluctuating air temperatures causing increased freezing and possibly reducing the infiltration capacity of the soil (Isard and Schaetzl, 1998; Iwata et al., 2011). How and where the snowpack melts is also an important component that has been overlooked primarily due to data limitations. Snow depth is heterogeneous at fine spatial scales, which is not possible to capture using point-scale station observations. Depth is also a poor indicator of the snow-water equivalence (SWE) of the snowpack, which can remain static as depth decreases and the snowpack goes through melt and refreeze cycles, causing the snow to increase in density. There is also a lack of available information on changes to groundwater recharge in warmer winters. Historical records of groundwater data beyond the last couple of decades are sparse and suffer from the same lack of spatial coverage as snow data. However, groundwater is a critical component of regional water budgets that should be examined to fully understand the changes to winter-spring hydrology (Sykes et al., 2016). Ford et al. (2020) estimated net recharge using a simple water balance equation and found that net recharge appears to increase in warmer winters in Michigan; however, the study was over short temporal and spatial scales that are insufficient to analyze climate-influenced trends. One solution to quantify the groundwater system’s 68 response to changes to melt, and to address limitations of this study, would be to develop process-based models of the hydrologic system, which can provide more detail on complex changes to the winter-spring hydrology. The changes to snow observed in this study are already affecting various facets of the region’s economy. Burakowski and Magnusson (2012) examined how climate change effects have impacted winter tourism economics in the United States and found that in years with lower snowfall from 1999-2010 most states with significant ski industries saw decreased visits, with a difference of $1.07 billion from high snowfall years. In the eastern US states, this economic difference ranged from -$51 million to -$40 million. This correlates with trends in other snow tourism related industries such as snowmobiling, which slightly declined from 2004-2010, with projected losses near $400 million annually in the Northeastern US by the end of the century if snow season declines continue under high emissions scenarios. Chin et al. (2018) reached similar conclusions about winter tourism in the Great Lakes region, with winter tourism businesses projected to decrease by 0.3-22% under high emissions forecast scenarios; this could lead winter tourism in many of these areas to become economically non-viable. Environmental and ecosystem responses to changes in snowmelt hydrology are also of concern. Increases in bare ground days and decreases in snow cover, especially in the more northern regions, could further exacerbate the warming climate as snow albedo decreases (Hayhoe et al., 2007). Decreases in snow cover and depth could have adverse effects on different vegetation species across the region as the reduction in snow could enhance freeze-thaw cycles leading to increase root damage and decrease frost protection (Perfect et al., 1987). Finally, the response of stream ecosystems to the shift in streamflow timing could negatively impact native aquatic species such as brook trout, which have significantly negative responses in spawning 69 population survival due to changes in spring stream temperatures and flow rates (Kanno et al., 2015). These ecological effects have been observed in the eastern US, but global implications are likely similar for other high-latitude regions. 5. Conclusions The results of this study confirm the stated research hypotheses that warmer winters are becoming more common, snowfall amounts and snowpack depths have decreased along with increased bare ground days, higher winter streamflows and earlier, reduced spring peak flows. Winters were classified as warm or cool at the basin scale using GHCN temperature data, then GHCN snow depth and snowfall data were examined along with USGS daily streamflow data. On average there were 3.57 more warm winters per decade in the 2010’s than in the 1960’s. These warmer winters had 50.1 cm less snowfall, a snow season length that was 24 days shorter, and an increase of 34 more bare ground days over this six-decade period. This correlates with changes to streamflow, including higher winter flows, three day earlier spring peak flows and 14 day earlier center of volume of discharge (CDV) times. The magnitude of difference in snow and hydrology between winter types was regionally-variable, with northern and western regions of the study showing the largest snowpack differences, northern and eastern regions having the largest differences in seasonal streamflow, and the eastern regions the largest difference in cumulative flows. These findings help fill an important gap in the literature regarding snowmelt hydrology in non-alpine settings and will provide a basis to improve projections of future changes to the high-latitude temperate snow seasons across the globe as winter temperatures continue to rise as a component of climate change. Despite extensive and relatively fine temporal resolution of weather station and stream gauge data used in this study, the spatial resolution is still coarse and the relationships correlative 70 in nature. Future studies should focus on causative analyses at higher resolution using process- based modeling of the hydrology related to snowmelt. Historical observations of streamflow and precipitation could be utilized in model calibration. Such a calibrated model could elucidate changes in hydrologic processes, which are not directly observable with existing data such as shallow groundwater recharge. Acknowledgments The text of this chapter is a reprint of the material as it appears in Science of the Total Environment in a manuscript first published June 18, 2021 (doi:10.1016/j.scitotenv.2021.148483) coauthored by Anthony Kendall and David Hyndman. This research was supported by USDA-NIFA grant 2015-68007-2313, NOAA grant NA12OAR4320071, and NASA grant NNX11AC72G along with Michigan State University’s Environmental Science and Public Policy Program and the Department of Earth and Environmental Sciences. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of USDA NIFA, NASA, or NOAA. 71 APPENDIX 72 Figure A.3.1. Additional study region map. HUC-2 basins used for grouping in this study, as well as the 1,751 USGS gages that provided streamflow data. 73 Figure A.3.2. Regional winter type precipitation. Median monthly precipitation for HUC-2 basins in warm and cool years. Generally, precipitation amounts are low in the winter months (DJF) and increase to their maximum amounts in the summer (JJA). Slight differences exist in precipitation amounts between different year types, but most basins show similar patterns regardless of year type. 74 Figure A.3.3. Regional winter temperatures. Median winter (DJF) temperatures for HUC-4 basins in the study region. Median temperatures range from -13.8°C to 4.8°C with temperatures increasing southward along a latitudinal gradient. 75 Figure A.3.4. Methodology flowchart. Flowchart visualizing the methodology applied in this study. 76 Figure A.3.5. Regional metric score slopes. Map showing the slope of the metric score per snow season for each HUC-4 in the study area. Only 3 of the basins had a negative slope, with a minimum of -0.0025. Most basins had an increasing metric score through time, especially along the Mid-Atlantic coast. 77 Figure A.3.6. Regional winter vs annual discharge. Each panel is a scatterplot of the mean annual basin yield in winter months (DJF) divided by the mean basin yield in all other months and then converted to a percent with a trendline fitted to the data. The basin names are: a) Souris- Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic. All basins except the two southernmost basins (h. and i.) have increasing trends through time. 78 Figure A.3.7. Regional basin yield difference. The difference in warm and cool winter cumulative basin yields. Median cumulative basin yields were extracted across winter type and HUC-2 basins (Fig. 7), then the cool winter cumulative basin yield was subtracted from the warm. Panels are approximately arranged by latitude of the HUC-2 for the: a) Souris-Red-Rainy, b) Great Lakes, c) New England, d) Upper Mississippi, e) Ohio, f) Mid-Atlantic, g) Missouri, h) Arkansas-White-Red, and i) South Atlantic Basins. A small inset of the HUC-2’s with that panel’s basin filled in is included for spatial referencing of the data. Note the different y-axis scales for each row. All but two of the basins (g. and h.) have higher cumulative basin yields in warm winters, with the increase starting during the spring snow melt period then starting to decrease at the tail end of the water year as streamflow becomes more baseflow dominant. 79 Winter Type Grouping Variable Warm Cool Metric Score 1.06 -0.93 Winter Temp (°C) -1.5 -5.5 Climate Annual Precipitation 92.6 88.7 (cm) Seasonal Snowfall (cm) 80.0 130.1 Max Snow Depth (cm) 23.5 37.9 Max Depth (DOWY) 111 123 Snowpack Season Length (days) 86 110 Bare Ground Days 202 168 Total Basin Yield (cm) 24.1 21.2 Max Basin Yield (cm) 0.12 0.14 Streamflow Max Basin Yield 192 195 (DOWY) CDV (DOWY) 172 186 Table A.3.1. Winter type statistics. Climate, snow and streamflow statistics generated across the entire study area for both warm and cool winters. All statistics were calculated as the median across winter types from annual HUC-4 statistics. Warm winters have more annual precipitation but less of that precipitation as snow. Snow in warm years is shallower, persists less and melts earlier leading to lower and earlier peak streamflow. 80 Ranked KS Grouping Variable Sum Annual <0.001 <0.001 Precipitation Peak Snowfall <0.001 <0.001 Climate Peak Snowfall Day of Water <0.001 0.009 Year Total Snowfall <0.001 <0.001 Peak Snow Depth <0.001 <0.001 Peak Snow Depth Day of Water <0.001 <0.001 Snowpack Year Season Length <0.001 <0.001 Bare Ground Days <0.001 <0.001 Seasonal Basin 0.051 0.089 Yield Peak Basin Yield <0.001 <0.001 Streamflow Peak Basin Yield Day of Water 0.007 <0.001 Year CDV <0.001 <0.001 Table A.3.2. Statistical test results. Results of the statistical tests of difference applied to warm and cool winter datasets. Only the seasonal basin yield totals were not determined statistically different at a 95% confidence threshold. All other datasets examined were statistically different according to both tests applied. The two-sided Student’s T-test was applied to the peak snowfall day of water year data (p-value < 0.001), as that was the only dataset with a normal distribution according to the Shapiro-Wilkes test. 81 CHAPTER 4: SNOWMELT DYNAMICS FROM WARMER WINTERS INCREASES GROUNDWATER RECHARGE IN MICHIGAN Abstract Earlier snowmelt and decreased snow amounts resulting from warmer winter temperatures with more bare ground days is well documented in many regions of the globe that receive significant annual snow amounts. The associated hydrology changes have significant ecological, infrastructure and economic effects. We examine changes to snowmelt hydrology between 2000 and 2020 across Michigan using a fully-distributed process-based landscape hydrology model to quantify changes to snowmelt and shallow groundwater recharge in warmer winters. We found that warmer winters in Michigan experience decreased snow amounts and earlier melt timing compared to cooler winters, with an average recharge of 3.6 cm less water in snow in warm versus cool winters and there is significantly less time that snowpack covers the surface. Differences in recharge between winter types were regional. In all regions of the state, recharge in warm winters was higher in winter months and lower in spring months compared to cool winters. In warm winters relative to cool winters, the Lower Peninsula had 3.4 cm more recharge, while the Upper Peninsula had 2.4 cm less recharge. The difference in recharge patterns over the course of the water year display greater differences in the northern regions than in the southern, suggesting the amount of overall snow received in a region plays a role in the how strongly recharge is affected by warming temperatures. 1. Introduction As climate change continues to warm the globe, winter temperatures in many mid-to-high latitude areas are increasing leading to changes to winter precipitation and snowmelt (Adam et al., 2009; Hyndman, 2014). Until recently much of the literature on these winter changes have 82 focused on alpine or arctic settings, leading to a knowledge gap in winter hydrology changes in non-alpine, seasonal snow cover environments. Although there is some emerging research on this topic, there is a strong need for additional research given the importance of snow regions such as the Great Lakes Basin, with its large surface water and groundwater reservoirs. The changes to winter hydrology in these regions can lead to significant ecologic and socioeconomic challenges as winter snowpacks decrease, spring streamflow patterns shift, and groundwater recharge to water table aquifers is altered. Winter precipitation is shifting to more rain and a smaller ratio of snow to total precipitation (S/P ratio) (Berghuijs et al., 2014; Feng and Hu, 2007; Ford et al., 2020; Hayhoe et al., 2010; Hodgkins et al., 2007; Huntington et al., 2004; Javed et al., 2019; Mankin and Diffenbaugh, 2014; Mote, 2003). These increased temperatures have led to decreased snowfall amounts, snowpack thickness, and snow cover, as well as increased numbers of bare ground days and earlier melting of the snowpack (Burakowski et al., 2008; Campbell et al., 2011; Clow, 2010; Dyer and Mote, 2006; Ford et al., 2020; 2021; Gan et al., 2013; Hamlet et al., 2005; Hayhoe et al., 2007; Hodgkins and Dudley, 2006a; Hodgkins et al., 2007; Ishida et al., 2019; Jefferson et al., 2008; Johnson and Stefan, 2006; McCabe and Wolock, 2010; Mote, 2003; Stewart et al., 2004; Suriano et al., 2019). The changes to melt timing and amounts have effects that propagate through the rest of the hydrologic system. Winter and spring streamflow in these snowy regions are changing in response to these melt changes. Specifically earlier and lower melt amounts have led to higher winter flows and lower spring flows, earlier peak flows, reduced peak flow amounts (and earlier center of volume (CV; also referred to as the center of discharge volume or “CDV”) of streamflows which is defined as the timing of the 50th quantile of flow in that water year (Campbell et al., 2011; Clow, 83 2010; Ford et al., 2020; 2021; Hayhoe et al., 2007; Hidalgo et al., 2009; Hodgkins et al., 2003; Hodgkins and Dudley, 2006b; Jefferson et al., 2008; Stewart et al., 2004). These hydrologic changes may not be limited to surface flows, with potentially increased recharge of shallow groundwater due to periodic snowpack melting leading to increased infiltration of the subsurface by melt (Ford et al., 2020). However, the research into groundwater changes in response to changing melt has been limited. This research addresses some of the knowledge gaps in snowmelt hydrology research. Specifically, our objectives are to simulate snow and hydrologic changes in the state of Michigan, which is a non-alpine setting in the Great Lakes Region of the United States that experiences variable winter snowfall amounts. Using a simulation model, we then quantify the changes in shallow groundwater recharge in warm vs. cool winters. Following up on earlier research in Ford et al. (2020), we hypothesize these warmer winters will correlate with increased recharge amounts when compared to cool winters with larger snow amounts and a more “typical” melt regime of one large melt in the spring compared to numerous periodic melts throughout the snow season. We will first describe the methodology implemented to test these hypotheses by describing the study area and period, along with the data sources and simulation model used to examine these hydrologic changes. Results from the model will then be presented and compared to previous research findings. The results are summarized, and future research questions are presented. 84 2. Materials and Methods 2.1 Study Region The state of Michigan in the center of the Upper Great Lakes Basin in the northeastern United States and encompasses a diverse landscape, geology, and climate (Fig. 4.1). The two peninsulas of the state are distinct in their land cover and geomorphological characteristics. The Lower Peninsula is dominated by large tracts of agricultural land interspersed with sprawling urban centers. These anthropogenically managed land cover types diminish to the North in the Lower Peninsula, giving way to primarily forests. Nearly all the land in the Lower Peninsula is underlain by deep glacial outwash deposits, which are composed primarily of sands and silts covering the various sedimentary rocks that make up the Michigan Basin. This contrasts with the Upper Peninsula, which is almost entirely remote forested wilderness, with very few urban areas or agricultural fields. The Upper Peninsula’s geology consists of very thin soils underlain by low hydraulic conductivity Precambrian bedrock. 85 Figure 4.1. Study region with the model boundary outlined with a thick black line, GHCN snow stations used for the “sliced” model cells, and model snowpack thickness assessment, HUC-10 basins used for spatial aggregation, and the larger colored regional designations used for visualization. The snowfall patterns also vary between and within the two peninsulas. Generally, the Upper Peninsula receives far more annual snowfall than the Lower Peninsula; some areas receive up to a half meter of snow water equivalent (SWE) totals (Ford et al., 2020). The Lower Peninsula has a latitudinal gradient in snowfall amounts, with the northern regions receiving more snow than the southern. There is also an east-west gradient in the Lower Peninsula driven by the Lake Effect phenomena where evaporated water from Lake Michigan condenses shortly after landfall, enhancing snowfalls near the peninsula’s western shore (Andresen, 2012). This results in the southeastern portion of the Lower Peninsula receiving a flux of less than 20 cm of 86 SWE and the northwestern regions receiving a flux more than 40 cm in some years (Ford et al., 2020). The Upper Peninsula also experiences Lake Effect snow generated by moisture from Lake Superior, causing the northern regions to receive more snow. The hydrology also varies across the state, driven by differences in landscape, climate, and geology. The state is dissected by a plethora of streams that drain surface runoff to the Great Lakes with perennial flow sustained by substantial baseflow contributions. The major hydrologic differences across the region lie in the groundwater characteristics. The Lower Peninsula contains extensive confined and unconfined aquifers. Many of the sedimentary bedrock formations contain aquifer units, and the deep, porous soils lend themselves to a robust unconfined shallow groundwater system. The Upper Peninsula also has vast amounts of shallow groundwater, but the thinner soils and impermeable bedrock don’t allow as much groundwater storage per unit area as the Lower Peninsula. 2.2 Landscape Hydrology Model (LHM) The simulations for this study were conducted using the Landscape Hydrology Model (LHM) (e.g., Hyndman et al., 2007; Kendall, 2009). LHM is a fully discretized, process-based simulation model. It is fully modular, and for this study is coupled with MODFLOW-2005 for groundwater simulations (Harbaugh, 2005). This simulation runs at hourly timesteps with a 1-km grid and four hydrologic zones: surface, root zone, deep unsaturated zone and saturated zone. LHM calculates each portion of the full water balance using the equation: ∆𝑆 = 𝑃 − 𝑇 − 𝐸 − 𝑃𝑐 + 𝑇𝑟 − 𝐸𝑥 − 𝑅 where: ΔS is the change in surface soil moisture storage, P is the watershed available precipitation, 87 T is the transpiration, E is the evaporation, R is runoff, Pc is percolation beneath the root zone, Tr is throughflow, and Ex is outflow from the cell. The snow model of LHM is based on the Utah Energy Balance (UEB) snow model (Tarboton and Luce, 1997). UEB is an energy budget snow model with a single layer of snowpack that primarily outputs snowpack SWE. The state of the snowpack is governed by the energy balance equation: 𝑞𝑁𝑒𝑡 = 𝑞𝑅𝐴 − 𝑞𝑆𝑢𝑟𝑓 − 𝑞𝐺 + 𝑞𝑀𝐼 + 𝑞𝐴𝑑𝑣 Where: qNet is the net energy budget of the snowpack, qRA is the net radiation absorption from short and longwave radiation, qSurf is the conductive heat flux with the atmosphere, qG is the conductive heat flux with the substrate, qMl is the latent heat of fusion of melt and qAdv is the advective heat flux from incoming precipitation and melting of the surface layer. While UEB provides a foundational framework for snowpack simulation, we altered several aspects of the model for this study. The primary modification was an additional calculation for snowpack thickness in addition to SWE to compare simulated and observed snow data, since many stations only collect snow depth measurements, and few stations in this region have SWE data available. To accomplish this, the snowpack model was updated to include variables for the snowpack liquid holding capacity and the fraction of ice/water within the 88 snowpack. This allowed for melt at the surface of the snowpack to be retained within the pack until the liquid holding capacity is met, with excess water drained from the snowpack as melt. The snow thickness after melt for each cell is generated by the equation: 𝑆𝑡 = 𝑆𝑡−1 − 𝑆𝑚 ⁄𝜌𝑤 where: St is the snowpack thickness, St-1 is the previous timestep thickness, Sm is the melt, and ρw is the density of water. Snow thickness after accumulation is calculated using the equation: 𝜌𝑠 𝑆𝑡 = 𝑆𝑡−1 + 𝑆𝑤 ∗ 𝜌𝑤 where: St is the snow thickness, St-1 is the previous timestep thickness, Sw is the incoming snow in SWE, ρs is the density of the snowfall and ρw is the density of water. To increase processing speed and efficiency during model runs, a “slice” model grid was developed. Instead of running the simulation across the entire model domain, 285 cells were first selected for simulation that contain weather station sites with observational snow data (Fig. 4.1). This increases the simulation speed for the entire 20-year study period, cutting down the time to results from multiple days for the entire grid down to several hours. Model outputs from both the full domain grid and the slice grid primarily used in this study’s analyses are snowpack thickness, snowpack SWE, melt from the snowpack and shallow groundwater recharge. 89 2.3 Data Sources A variety of observations and model outputs were used for model inputs and calibration. The snow module of LHM was calibrated to two primary sources: empirical snow data from the Global Historical Climatology Network (GHCN) weather station data (Menne et al., 2012), and Daily observations of snow depth at 285 stations across Michigan retrieved using the rnoaa software. The location of these stations were used to identify the grid cells used in the sliced model (Chamberlain et al., 2021, Fig. 4.1). In addition to the observed snow data, outputs from the National Oceanic and Atmospheric Administration’s (NOAA) Snow Data Assimilation (SNODAS) model were compared to our model (NOHRSC, 2004). SNODAS provides daily SWE values across the study area on a 1-km grid back to 2003. SNODAS utilizes downscaled numerical weather model outputs, assimilated with observational snow and climate data from weather stations and satellite data (Barret, 2003). This model was chosen for snow calibration because, in addition to the coverage and spatial resolution, the data assimilation model adequately represents snowpack components such as SWE and snow depth in relatively flat, non- alpine settings such as Michigan (Clow et al., 2012; Hedrick et al., 2015). 90 Figure 4.2. Selected inputs for LHM include: a) mean water year SWE across the study region from NLDAS-2 forcing data (Xia et al., 2012); b) LULC classifications for the region simplified from NLCD data (USGS, 2011); c) the saturated hydraulic conductivity of sediments in the top soil layer from gSSURGO (Soil Survey Staff, 2020); and d) topography from the NED DEM (USGS, 2020). Inputs to LHM are divided into four categories depending on which parts of the model they influence and the nature of the data: climate, landscape/topography, static inputs and observational inputs. Each category has variable inputs, which change depending on the aim of a 91 particular simulation. Climate inputs such as air temperature and precipitation come from NLDAS-2 forcing data (Xia et al., 2012; Fig. 4.2a). Observational data came from two sources: streamflow data came from the US Geological Survey’s stream gaging network (USGS, 2018) and groundwater level data from the Michigan Department of Environmental Quality’s (MDEQ) “WelLogic” database, which contains static water levels for wells across the state at the time of well installation (MDEQ, 2008). For landscape inputs we used: a 30 m (1 arc second) digital elevation map (DEM) from the USGS National Elevations Dataset (NED) (USGS, 2020; Fig. 4.2d), leaf area index (LAI) values from MODIS (Myeni et al., 2015), and land use-land cover (LULC) data from the National Land Cover Dataset (NLCD) from 2001 to 2016 (USGS, 2016; Fig. 4.2b). The NLCD data was categorized using a modified land cover classification scheme of Anderson classes (Anderson et al., 1976). Static inputs include: gSSURGO soil data from the US Department of Agriculture (Soil Survey Staff, 2020; Fig. 4.2c), wetland and lake depth information from the National Wetlands Inventory (NWI) (US Fish & Wildlife Service, 2018), hydrography data from the National Hydrography Dataset (NHD) (USGS, 2019), initial bedrock hydraulic conductivity (K) and specific storage (Ss) estimates from the Lake Michigan Basin Model (Feinstein et al., 2010) plus K and Ss for surficial deposits based on the WelLogic dataset (MDEQ, 2008). 2.4 Model Calibration Several calibration steps were conducted to improve the accuracy of the snowpack simulation. The first step involved comparing the slice model grid daily values of snowpack thickness and SWE to station values of daily snow depth and SNODAS SWE values. SNODAS comparisons to LHM were cell-to-cell comparisons determined by the proximity of cell centers since the two model grids don’t exactly overlap. These values were primarily used to calibrate 92 the depth of ground heat exchange, which is a scalar value (denoted as ‘de’ in the UEB code). This value was originally set as 0.6 m in UEB, but was changed to 1.0 m for these simulations based on the calibration. There are assumptions about using another gridded simulation product such as SNODAS to calibrate a separate simulation such as LHM, namely the accuracy of SNODAS in representing the snowpack. SNODAS has been shown to be reasonably accurate in areas such as Michigan where the landscape is relatively flat, resulting in negligible influence from aspect and gravity effects on snowpack movement and melting (Clow et al., 2012; Hedrick et al., 2015). Observational data such as at GHCN stations is far more useful for improving the accuracy of the model, but the sparse spatial distribution of this point data limits its usefulness for a model of regional scales such as this, and many of the observational weather data may not be representative of the area as they are typically obtained in flat areas without any canopy cover. Thus, to improve the accuracy of the LHM simulation both data types were utilized in calibration of the snow model. 2.5 Analysis Methods Once the simulation was complete, several analyses were conducted to quantify changes to melt hydrology in warmer versus colder winters. First, full grid simulation outputs were aggregated to 10-digit Hydrologic Unit Code (HUC-10) basins (Seaber et al. 1987; Fig. 4.1). The basin aggregated values for each water year were then classified as “warm” or “cool” using a multimeric analysis developed where six temperature related metrics for each snow season (defined as Oct-May) were compared to the overall basin norm as described in Ford et al. (2021) (Fig. A.4.1). Using these winter type classifications, warm and cool winter simulations were compared to observe differences in annual basin aggregations of SWE, depth, melt and recharge. Each winter type dataset was tested for normality using the Shapiro-Wilk Normality Test 93 (Royston, 1995), and datasets confirmed as non-normal distributions were analyzed using the Mann-Whitney (Wilcoxon) ranked sum test (Hollander and Wolfe, 1973) and the two-sample Kolmogorov-Smirnov (KS) test (Conover, 1971). The p-values for these statistical tests were then evaluated with a 95% confidence threshold. For better visualization, results were aggregated across six regions of the state: the Upper Peninsula was bisected latitudinally, and the Lower Peninsula was broken into four quadrants (Fig. 4.1) based on differences in annual snowfall amounts. All postprocessing of the simulation results and subsequent analyses and visualizations were performed in the “R” (R Core Team, 2021) programming software. 3. Results 3.1 Snow Model Performance The simulation of the snowpack by LHM performed better in terms of snowpack SWE than depth. When aggregated across the entire study area LHM’s SWE was lower than SNODAS by several centimeters, especially in the late winter-early spring months (Fig. 4.3a.). However, this varies regionally, as the difference between LHM-SNODAS typically varies by less than 1 cm for all but the northernmost region in the study area (Fig. 4.3c). 94 Figure 4.3. LHM snowpack model compared to SNODAS and GHCN data. All data is aggregated across the entire study region and all years as a) mean monthly SWE averaged across the study region; b) mean monthly depth of model slice cells with LHM overpredicting compared to station data; c) the difference in LHM-SNODAS monthly SWE-all regions except the Northern UP have 1 cm or less difference; d) the LHM-GHCN difference in depth showing the northern regions having a smaller difference than the southern regions. In general, LHM overpredicts the snowpack depth in all but the first couple months of winter by tens of centimeters. (Fig. 4.3b). In contrast to the SWE simulations, the difference between LHM depth and GHCN depth is most significant in the Lower Peninsula regions, with the Southeastern Lower Peninsula reaching a difference of 44 cm in March, while the two Upper Peninsula regions have less difference in the opposite direction (Fig. 4.3d). 95 3.2 Snowpack Analysis Figure 4.4. HUC-10 SWE totals for each complete water year in the study period. SWE increases with latitude in nearly all water years. Some years (2001, 2014, 2019) are colder in all regions, while other years show more regional variation. Snowpack simulations capture regional gradients as well as temporal differences. In all winters within the study period, the largest annual snowpack SWE amounts occur in the Northern Upper Peninsula basins, with these basins reaching over 30 cm of SWE in the coldest years (Fig. 4.4) and having annual maximum depth averages over a meter in cool winters (Table 4.1). By contrast, the Lower Peninsula has much lower snow amounts, especially in the southern regions where annual SWE totals are typically ten centimeters or less, and annual maximum depths as low as 14 cm in warm winters (Fig. 4.4; Table 4.1). Melt values show similar patterns, with the highest annual melts in the northern regions and during cooler winters. 96 Region N UP S UP NW LP NE LP SW LP SE LP Winter Type Warm Cool Warm Cool Warm Cool Warm Cool Warm Cool Warm Cool Total Precipitation 85.9 85.1 84.5 82.2 87.0 85.8 81.9 81.2 93.0 92.9 89.8 86.3 (cm) Max Snow 83.2 132.8 65.7 109.0 59.8 90.3 43.0 68.5 26.5 67.5 21.5 58.2 Depth (cm) Max SWE 10.6 14.2 9.0 12.6 8.9 10.3 6.7 8.2 4.5 9.5 3.6 8.3 (cm) Total Melt 17.3 21.8 14.1 18.7 15.8 17.4 14.3 16.6 10.3 15.3 9.9 13.3 (cm) Total Recharge 26.3 29.7 23.2 24.6 41.8 39.3 29.7 28.3 26.9 23.2 21.8 16.0 (cm) Table 4.1. Regional winter type statistics for selected hydrologic model outputs. Precipitation totals are similar between winter types, while snow depth and SWE vary significantly in all regions. Recharge totals are higher in warm winters in southern regions while lower to the north. Figure 4.5. Mean monthly SWE in warm and cool winters for a) the northern UP, b) the southern UP, c) the northwestern LP, d) the northeastern LP, e) the southwestern LP and f) the southeastern LP. In all regions cool winter SWE peaks higher and later in the season. The difference between cool and warm winter snowpacks is not just in the annual totals, but also in the evolution of the snowpack throughout the winter. In all regions, snowpack SWE 97 peaks earlier in warmer winters, and the peak is a lower (Fig. 4.5). The difference between winter types is smallest early in the winter, especially in the northern Lower Peninsula regions, and becomes larger as the snow season progresses, only decreasing in the late spring when the snowpacks become greatly diminished. In the Lower Peninsula regions, the snowpack in cool winters persists a month longer than in warmer winters, where the snowpack has usually completely melted by the end of March. 3.3 Recharge Analyses Figure 4.6. Mean simulated recharge by winter type in: a) the northern UP, b) the southern UP, c) the northwestern LP, d) the northeastern LP, e) the southwestern LP and f) the southeastern LP. Winter month recharge is higher in warm winters in all regions, but peak recharge in cool winters is higher in all but the southeastern LP. The groundwater recharge time series for each region shows distinctive differences between warm and cool winters. In all regions, winter recharge is higher in warm winters, with the largest difference in the southern regions (Fig. 4.6). Cool winter recharge is higher in the spring, and then remains higher for the rest of the water year, with the difference between winter 98 types diminishing later in the water year. These late season differences are most prominent in the UP regions, and smaller in the two southernmost regions. Recharge peaks 1-2 months earlier in warm winters in the Lower Peninsula regions. The timing in recharge peaks are the same in the UP regions, but the cool winter recharge peaks are several centimeters higher, while the recharge peaks in the LP regions are similar. Cumulative recharge better shows the difference between winter types by region, with warm winter recharge remaining higher than cool winter recharge throughout the entire water year in the LP, while the UP has very similar recharge amounts in both winter types until the late spring and summer months, when cool winter recharge increases relative to warm winter recharge (Fig. 4.7). Figure 4.7. Cumulative monthly warm and cool winter recharge for a) the northern UP, b) the southern UP, c) the northwestern LP, d) the northeastern LP, e) the southwestern LP, and f) the southeastern LP. In the UP regions, cool winter recharge totals are higher, but warm winter total recharge is higher in LP regions, with larger differences moving southward. 99 3.4 Statistical Analyses When applied to the HUC-10 basin annual averages of the different hydrologic variables examined, the Shapiro-Wilk test confirmed all datasets as non-normal. Evaluating the statistical tests with a 95% confidence threshold, only one variable did not meet this confidence threshold for difference between the datasets was runoff. The p-value from the Mann-Whitney Ranked Sum test for runoff was 0.29 (Table 4.2). Using the Kolmogorov-Smirnov analysis the runoff datasets did meet the confidence threshold with a p-value of 0.03. All other datasets evaluated resulted in p-values below 1% for both tests applied. Ranked Variable Sum KS Precipitation <0.001 <0.001 SWE <0.001 <0.001 Depth <0.001 <0.001 Melt <0.001 <0.001 Recharge <0.001 <0.001 Table 4.2. P-values for statistical tests of difference for selected output variables. All tests met the 95% confidence threshold except for runoff in the Mann-Whitney test, which is italicized. 4. Discussion Recharge flux is higher in the Lower Peninsula, whereas in the Upper Peninsula, recharge flux is greater in cool winters. These regional differences are likely heavily influenced by differences in Land Use-Land Cover and latitudinal differences in snow amounts. Urbanization increases from north to south, and in conjunction with the lower snow amounts, relative to the other regions, likely leads to similar recharge in both winter types as snow amounts are closer and there is less permeable soil to infiltrate. In the Upper Peninsula regions, the heavy snow amounts, especially in cool winters, likely lead to higher recharge in cool years because there is significantly more water introduced into the system. Even in warm winters these regions see as much snow as the southern regions do in cool winters (Fig. 4.5). 100 Despite the difference in regional results, the changing snowpack dynamics in warmer winters leading to increased recharge can be inferred based on the timing of when that recharge occurs. In the northern portion of the Lower Peninsula regions, warm winter recharge peaks in March compared to April in cool winters, which is consistent with the earlier melt timing (Fig. 5.6). These peaks are also lower in warm winters. In the southern regions of the Lower Peninsula, the recharge peaks are similar in both winter types, as is timing, but the amount of recharge flux in the late fall and winter months is so much larger than in cool winters that the difference in recharge flux between warm and cool winters is the highest in the study region. These recharge patterns are very similar to the net groundwater recharge estimates found in Ford et al. (2020), which used a simple water balance equation and found net recharge was higher in warm winters in the south, but higher in cool winters in the Upper Peninsula. As anthropogenic climate change leads to increased frequency and intensity of warmer winters (Ford et al., 2021), these melt hydrology changes could have significant infrastructure and ecological effects. Higher winter melt amounts and recharge in areas with shallow water table levels could lead to increased incidences of winter and early spring flooding. The increased recharge in warmer winters could lead to changes in spring streamflow which is a critical spawning time for many freshwater fish species. Both changes to recharge could also affect local and regional water supplies in Michigan, as many localities rely on groundwater reservoirs for domestic use. The results of this study have increased confidence in the model utilized to observe these snow changes when compared to previous studies that rely entirely on observational data and reach similar conclusions. Despite this, as with any simulation model, a number of uncertainties exist. While the basin scale is necessary to visualize and interpret results over such a large area, 101 snow is extremely heterogeneous at fine spatial scales, leading to potential inaccuracies at finer resolutions. The model also only simulates the snowpack as a single layer, which is not representative of snowpack evolution where melting and refreezing within the snowpack followed by additional accumulation leads to distinct layering within the snowpack. This likely contributes to the difference in modeled snowpack thickness relative to the observations. Finally, a model is only as good as its inputs. While the inputs utilized in this study are well established and recognized as reasonably accurate, there are differences in spatial and temporal resolutions compared that may have led to some errors. Future studies will require further refinement of the snow model in LHM. The thickness simulation will need further calibration to better match observed depths. Including multiple layers within the snowpack would greatly increase the snow model’s accuracy, both in depth and simulation of melt from the snowpack. The regional differences observed in this study also deserves further examination to see if the proposed interpretations related to LULC and snow amounts are correct. 5. Conclusions Using the Landscape Hydrology Model, this study simulated snow, recharge, and runoff across the state of Michigan from 2000-2020. The outputs from this simulation were categorized at the basin level into warm and cool winters based on temperature norms, finding changes to seasonal snowpack evolution and melt hydrology in warmer winters. Specifically, warmer winters led to decreased snowpacks which melted earlier and more frequently. These melt changes resulted in increased shallow groundwater recharge in the winter months. 102 Acknowledgments This research was supported by USDA-NIFA grant 2015-68007-2313, NOAA grant NA12OAR4320071, and NASA grant NNX11AC72G along with Michigan State University’s Environmental Science and Public Policy Program and the Department of Earth and Environmental Sciences. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of USDA NIFA, NASA, or NOAA. 103 APPENDIX 104 Figure A.4.1. Winter type classification flowchart. This flowchart visualizes the process for classifying winters as warm or cool and how that classification is applied to analyses results. 105 CHAPTER 5: A WARMER FUTURE: CLIMATE CHALLENGES FOR MICHIGAN’S SKI INDUSTRY Abstract Warming winter temperatures in mid-to-high latitude regions is changing the amount and timing of snowfall in these regions, with warmer winters leading to decreased snowfall, lower snowpack thicknesses, less persistent snow cover throughout the winter, and earlier snowpack melting. These changes could have major economic effects on the winter ski industry, which generates billions of dollars annually across the United States. In Michigan, the ski industry generated 10,889 jobs and $638.3 million revenue in 2009-2010 (Burakowski and Magnusson, 2012). This study examines the attendance difference between warm and cool winters in the Michigan ski industry. These snow differences provide a lens to examine climate responses. Across the midwest, ski slope visits decreased by an average of 387,000 visitors per year in warm winters compared to cool. At Shanty Creek in the northern portion of the Lower Peninsula of Michigan, water usage for snowmaking was 37,470 m3 higher in cool winters compared to warm and had 37 additional snowmaking hours due to the longer cool periods. 1. Introduction North America’s snowy regions are observing reductions in winter snow cover as winter temperatures increase (Burakowski et al., 2008; Chin et al., 2018; Dyer and Mote, 2006; Feng and Hu, 2007; Ford et al., 2020; Ford et al., 2021; Gan et al., 2013; Hamlet et al., 2005; Hodgkins and Dudley, 2006; Huntington et al., 2004; Ishida et al., 2019; Mankin and Diffenbaugh, 2015; Mote, 2003; Peacock, 2012; Suriano et al., 2019). These increased winter temperatures have led to declining snowpack thickness, earlier melting of the snow pack, more 106 bare ground days in winter and a shift to more precipitation falling as rain rather than snow (Berghuijs et al., 2014; Burakowski et al., 2008; Campbell et al., 2011; Clow, 2010; Dyer and Mote, 2006; Ford et al., 2020; 2021; Gan et al., 2013; Hamlet et al., 2005; Hodgkins and Dudley, 2006; Hayhoe et al., 2007; Huntington et al., 2004; Ishida et al., 2019; Javed et al., 2019; Jefferson et al., 2008; McCabe and Wolock, 2010; Mote, 2003; Stewart et al., 2004). These changes have far reaching effects globally, as 1.2 billion people rely on snowmelt to meet domestic and agricultural water consumption, and projected snowpack reductions resulting from a warming climate could lead to the loss of trillions of dollars (Sturm et al., 2017). These changes to snow amounts and timing have the potential for significant implications for the economies of regions that generate large amounts of revenue from winter recreation and tourism, in addition to the numerous hydrologic and ecological effects. There are over 450 ski slopes in the United States, down from over 700 in the 1980’s (Hagenstad et al., 2018). In the 2009-2010 winter season, $10.7 billion was generated for the US economy from skiing and snowboarding, with an additional $1.5 billion from snowmobiling (Burakowski and Magnusson, 2012). In 2015-2016 the amount generated by winter recreation in the US had increased to $20.3 billion (Hagenstad et al., 2018). These sales come from slope passes and rentals, as well as the associated tourism spending on lodging, food, and travel. Just on land owned by the US National Forest Service, where a quarter of the 470 ski slopes in the US are located, there are 2 million downhill skiers annually, making it the second most popular outdoor recreation on National Forest land after hiking; this adds $26 million annually to the US treasury (second only to timber) (Chapagain et al., 2018). Across North America there was an annual average of 72.6 million ski visits from 2011- 2016, a decline of 5.34 million from the previous five-year average (Knowles, 2019). There is a 107 high positive correlation between ski slope visits and snow cover amount, with the decreases in visits during low snow winters more significant than the increases in high snow winters (Hagenstad et al., 2018). From 2001-2016 there were 3.8 million more slope visits in high snow winters, compared to 5.5 million less in low snow winters (Hagenstad et al., 2018). This translates to an extra $692.9 million and 11,800 jobs in these high snow winters and a loss of $1 billion and 17,400 jobs in low snow winters (Hagenstad et al., 2018). From 1999-2010 lower annual snowfall amounts correlated with less slope visits in almost all states with ski areas (Burakowski and Magnusson, 2012). In Michigan, each inch of snow depth corresponded with an increase in daily lift ticket sales of 7-9%, with temperature having a negative correlation on sales (Shih et al., 2009). While more snow means higher ski demand, the influence of temperature is also important, with seasons with higher mean temperatures correlating with decreased ski trips (Chapagain et al., 2018). Sometime just the perception of poor ski conditions leads to declines in visits, with Hamilton et al. (2018) finding that ski visits in New England were correlated more heavily with snowfall amounts in nearby urban areas than the actual slopes. In addition to the snow pack and temperature influencing ski slope visits and therefore economic effects, the timing of when that snow occurs is also significant. Ski areas require at least 30 cm of base snow depth to function with at least 100 continuous days of at least 50% of slopes open to be profitable, with the 2 weeks around the Christmas holiday season being the most vital (Knowles, 2019). The average length of the ski season using these parameters peaked in North America in 2010, with the combination of later opening dates, earlier closing dates and increased mid-season closures contributing to the first decline in season length in over 30 years (Knowles, 2019). The later start dates are beginning to encroach on the holiday timeframe critical for ski revenue. A study by Wilson et al. (2018) using 50 years of data from the Hubbard 108 Brook Experimental Forest in New Hampshire found significant warming of winter temperatures during that period, with the most significant warming occurring between December 1-25. Skiing isn’t the only winter recreation suffering the economic effects of these reduced season lengths. In Vermont the amount of snowmobile registrations is declining, with regions in the Northeastern US potentially losing more than half the snowmobile season by the end of the century, likely making the sport unviable in many areas (Perry et al., 2018). The outlook for winter recreation opportunities and viability in North America is not encouraging given climate projections. Under different emissions scenarios, combined with projections in population growth, downhill skiing could lose 12-20% of current visits by mid- century and up to a 45% decrease by 2090 as the available area decreases and the slopes that remain become more crowded (Wobus et al., 2017). This same study examining 247 ski areas across the continental US found that a decrease in season length at all locations under different emissions scenarios, with more than 50% decrease in season length by 2050 and an 80% decrease by 2090 for some slopes. This would likely lead to a loss of ~$2 billion annually under high emission scenarios due to a roughly linear correlation between season length and ski visits (Wobus et al., 2017). The most significant reductions in the Wobus et al. (2017) study occurred in the upper Midwest and Northeastern US, with the smallest changes occurring in the alpine regions of the Western US like the Rockies and Sierras. These lower elevation slopes in the Eastern US could potentially see a near complete loss of current winter recreation activities by 2090 under high emissions scenarios. Nationally this corresponds with a decrease of downhill ski slope visits from 56 million annually to 19.8 million by 2090 under high emissions (Wobus et al., 2017). These outcomes are largely dependent on the greenhouse gas emissions scenario, with relatively 109 minimal losses to season length and skiable areas under low emissions, but under high emissions only 29 slopes in Quebec and the Northeastern US are open for at least 100 days and during the critical holiday season (Scott et al., 2021). There are similar findings for the Great Lakes region, where season length decreases 25-38% by midcentury depending on emissions, with much of that occurring at the beginning or end of the season and the holiday window reduced by up to 66% and snow depths during the holidays reduced by 50% or more by the end of the century under high emissions (Chin et al., 2018; Scott et al., 2021). This could lead to only 7-8% of ski areas in the Midwest to be financially viable by 2100 under high emissions. The ski industry is already adapting to the changes to the winter season, but ultimately it may not be enough. Artificial snowmaking technology and processes have advanced in recent decades, but it’s unclear if this will be able to keep up with projected changes. Currently 100% of ski areas in the Eastern US, 94% in the Midwest and 91% in the Rocky Mountain region utilize snowmaking (Ooi, 2017). Improvements made to snowmaking include high efficiency snowguns that have a 1:1 air/water ratio, gravity-fed snowmaking systems, upgraded computer automation and monitoring, more efficient air compressors, water coolers for lower ambient temperatures of the water used, flowmeters on snowmaking equipment for monitoring water use, and the use of reclaimed water (Ooi, 2017). However, snowmaking is limited both by the available water and energy necessary for operation and the winter air temperatures (Wilson et al., 2018). The source of water is a hard limit as only so much surface water can be utilized, and many of these streams may face decreased late winter/early spring flows with decreased snowmelt (Ford, et al., 2021; Wilson et al., 2018). One hectare of artificial snow at a depth of 30 cm requires 5,000-27,000 kWh of electricity and 600,000-1,500,000 L of freshwater (Scott et al., 2021). As temperatures continue to rise, reducing snow cover the demand for artificial snow will 110 increase; this will further stress these freshwater resources (Wilson et al., 2018). This will compound the temperature limitations of generating artificial snow, which needs a daily average temperature of -2 °C or less for adequate start-up of the snowguns (Wilson et al., 2018). Between the combined resource limitations and future temperature increases it is unlikely artificial snowmaking will remain a financially viable option under future high emissions scenarios (Knowles, 2019). The national and regional economic contributions of the ski industry, along with the increasing pressures resulting from climate change, necessitate a thorough understanding of the potential impacts of these changes for successful adaptation and mitigation methods. This study addresses this need by focusing on the winter climate and ski industry in the state of Michigan, which contains the most ski slopes in the Midwest, and the second most in the country (Scott et al., 2021; Shih et al., 2009). The specific objectives of this study are to examine changes to winter temperature and snow in Michigan and compare those to changes at ski slopes at a local and regional scale. First, the winter climate of Michigan in recent decades will be described to help contextualize these changes to the snowpack, and then the state of ski tourism in Michigan will be presented. These two components are then compared directly, and then the results utilized to discuss potential future effects. 111 Figure 5.1. Study region map. Average annual winter (DJF) snow depth for HUC-8 basins in Michigan from 2004-2020. Analyses in this research used the basins containing the ski areas shown. Shanty Creek Resorts is denoted in red. High snow depth areas in the state are in the Upper Peninsula and northern and western portions of the Lower Peninsula. 2. Methods 2.1 Study Region The state of Michigan is in the center of the Great Lakes Basin in the upper Midwestern United States (Figure 5.1). The two peninsulas of the state are surrounded by Lakes Michigan, Huron and Superior, enhancing precipitation across the state, particularly the western Lower 112 Peninsula and the northern Upper Peninsula. This enhanced precipitation is generated by the “Lake Effect” phenomena where winds across the lakes increase evaporation, followed by condensation of this moisture-laden air shortly after moving over the land surface (Andresen, 2012). This Lake Effect, along with the colder winter temperatures of the mid-to-high latitude state and it’s glacially influenced hilly terrane make Michigan a premier winter recreation destination in the Midwest. There are almost 100 ski areas located in the Midwest, with 40 of those located within Michigan, the most of any state in the region (Scott et al., 2021). During the 2015-2016 season, the Midwestern ski market generated $1.65 billion and approximately 27,000 jobs constituting around 11% of the total US ski market visits (Scott et al., 2021). In the winter of 2004-2005, Michigan had the second highest number of ski slopes in the country, attracting 40% of the states ski visitors from out of state (Shih et al., 2009). During the 2009-2010 winter, skiing and snowboarding created $3.5 billion in the Great Lakes region, and in Michigan $638.3 million was added to the state’s economy along with 10,889 associated jobs (Burakowski and Magnusson, 2012). The outdoor recreation and tourism industry is the third largest industry in the state of Michigan (behind auto manufacturing and agriculture), with skiing the most popular winter recreational activity in the state (Shih et al., 2009). Like many ski areas across the country, Michigan’s ski industry is likely to suffer negative impacts under a warmer future climate, however it is more resilient than more southern states in the region (Scott et al., 2021). Currently 76% of the ski areas in the Great Lakes can be considered economically viable, but that decreases to just 7-8% by the end of the century under high emissions scenarios (Scott et al., 2021). By the 2080’s the days with sufficient snow depths for winter recreation in the region could be up to a month shorter, with less than one month per 113 year with temperatures suitable for snow making (Chin et al., 2018). The financially critical holiday period is stable until the end of century under high emissions, where it’s reduced by up to 66% with a 50% or more decrease in snow depths during those weeks (Chin et al., 2018; Scott et al., 2021). While snowmaking capabilities may be able to counteract these changes in the early part of the century, by the end of the century under high emissions Michigan slopes’ snowmaking requirements increase 516% and the states skiable terrane is cut by more than half under high emissions compared to low (Scott et al., 2021). Considering the current importance of the ski industry to Michigan and the Great Lakes, such changes are likely to have significantly negative economic effects on the region. 2.2 Data Sources Several data sources of varying spatial and temporal scales were utilized in the analyses performed in this study. Temperature data used to classify the types of winters came from daily weather station observations in the Global Historical Climatology Network (GHCN) (Menne et al., 2012). Snow data was extracted from the National Oceanic and Atmospheric Administration’s (NOAA) Snow Data Assimilation (SNODAS) model, which is a daily simulation of snow depth (and other snowpack parameters) at 1-km gridscale (NOHRSC, 2004). These two data sources were spatially aggregated to Hydrologic Unit Code 8-digit (HUC-8) basins (Seaber et al., 1987). These watershed basins provide contiguous, non-overlapping basins of suitable size to capture regional trends. Only HUC-8’s containing Michigan ski areas were used in the analyses. The data were then examined during snow seasons, defined as October 1- May 31, from 2003-2020 (SNODAS products became available starting in 2003). This resulted in 21 basins for this study, primarily in the northern areas of the state. 114 In addition to the above data sources, further insight into ski slope trends came from annual statistics for the Shanty Creek Resort provided by Chief Executive Officer Pete Bigford. This data contains information on snowmaking time and resource usage for snowmaking from 2014-2021. Shanty Creek, located near Bellaire, Michigan in the northern Lower Peninsula has been in operation since the 1960’s with 53 ski runs (Shanty Creek, 2021). 2.3 Analysis Methods The analyses conducted here uses the year type classifications for winter described in Ford et al. (2020). GHCN temperature data were aggregated across the HUC-8 basins, then winters were classified as warm or cool relative to the norm across all years. Using those year type classifications snow data was examined for all basins containing ski slopes in Michigan. Specific snow variables examined in the year types include the mean and maximum snow depths, amount of bare ground days (defined as days with < 2.5 cm of snow depth) and the amount of days with minimum snow depth required for skiing (defined as a depth of 30 cm or greater). Special attention was given to the economically critical holiday period (12/22-1/2). In addition to being used for the year type classification, the aggregated temperature data was also examined for the mean snow season temperatures during the different year types. Because warmer winters will require increased snowmaking production, the number of days suitable for snowmaking was also quantified for the different year types, which is defined as the number of days with a mean daily temperature of -2°C or less (as described by Wilson et al., 2018). The conditions necessary for snowmaking are more affected by the wet bulb temperature rather than the commonly reported dry bulb temperature (Wobus et al., 2017; Pete Bigford, personal correspondence). Unfortunately, historical empirical data for the wet bulb temperature is not widely available, thus the dry bulb temperature must be used for a proxy in this study. 115 All of these variables were examined across all ski basins in the state, with detailed analysis of the basin containing Shanty Creek Resorts. Using data from the National Ski Areas Association (NSAA), annual regional ski area visits for the entire Midwest were also examined to determine the difference in ski area visits in warm and cool winters during the study period. The spatial temporal aggregation of the data and subsequent analyses performed in this study were all conducted in the “R” (R Core Team, 2021) programming software. 3. Results and Discussion The multimetric analysis of station temperature data across the state identified six warm winters and 3 cool winters during the 17-year study period. In all the basins containing ski areas in Michigan, warmer winters had significantly lower snow depths throughout the season when compared to cool winters (Figure 5.2). Mean winter depth across all basins was only 8.6 cm compared to 19.1 cm in cool winters, with a mean maximum snow depth that was 19.2 cm lower in warm winters (Table 5.1). The snow season was shorter in the ski basins in warm winters, starting later and ending 15 days earlier (the snow season here is defined as the first and last days of the season with a depth of 2.54 cm or greater). This resulted in only 2 days of 30 cm or greater depths in warm winters compared to 76 in cool, neither of which meeting the 100 days of 30 cm depth or greater for financial viability of ski slopes found in the Demiroglu et al. study (2016) (Figure 5.3a). For the economically critical holiday season, both warm and cool winters did not meet the 30 cm depth, with a mean depth of 12.4 cm in warm winters and 21.4 cm in cool winters for that 2-week period. 116 Figure 5.2. Winter type daily depth. Mean daily snow depth in warm and cool winters for all basins containing ski areas in Michigan (dashed line) and the Shanty Creek basin (solid line). The dashed horizontal line denotes the 30 cm depth necessary for adequate skiing, and the two vertical lines represent the economically critical holiday season. Water years start October first, which is Day 0 in the figure. Warm winters consistently have lower depths throughout the season and melt significantly earlier. 117 Region All Ski HUC-8's Shanty Ski HUC-8 Year Type Warm Cool Warm Cool Mean Depth (cm) 8.6 19.1 9 16.4 Max Depth (cm) 44.1 63.3 45.6 58.2 Bare Ground Days 152 107 147 101 Days with Minimum 25 70 25 60 Skiing Snow Depth Mean Temperature (°C) 3.4 -0.2 3.6 0.3 Days with Temperatures 66 99 61 91 Suitable for Snowmaking Average Visitors 6.671 7.058 NA NA (millions)* Water Use (m3) NA NA 158040 195511 Snowmaking Hours NA NA 710 747 Table 5.1. Winter type statistics. Warm and cool winter stats for all Michigan basins containing ski areas and the Shanty Creek basin. Warm winters on average have higher temperatures, lower snow depths, fewer days with the minimum depth required for skiing (30 cm), fewer days with adequate snowmaking conditions and less ski slope visitors. Shanty Creek’s basin’s snow shows similar warm and cool winter differences as the other basins of the state (Figure 5.2). Snow depths are lower in warm winters during the entire season, with a mean depth of 45.6 cm in warm winters and 58.2 cm in cool winters (Table 5.1). Mean maximum depth in warm winters was 45.6 cm, 12.6 cm lower than in cool winters. The snow season ended 26 days earlier in warm winters compared to cool, and only 6 days in warm winters met the 30 cm threshold, compared to 49 in cool winters. The holiday weeks saw a mean depth of 13.8 cm in warm winters, and only 14.7 cm in cool winters. The patterns of the snowpack match those found across all ski basins, with similar warm winter snowpacks. However, the differences between warm and cool winters are less pronounced than across the rest of the state. This is similar to the findings of Ford et al. (2020), where more southern basins in Michigan 118 showed less difference than the more northern basins, likely attributable to the overall less snow found in the south. There are significantly fewer days meeting the temperature threshold of -2 °C for snowmaking in warm winters compared to cool (Figure 5.3b). Across all basins, cool winters typically had 2-4 months of days at or below this temperature with a mean of 99 days (Table 5.1). In warm winters the amount of snowmaking days ranged from ~ 1 to 3 months, with a mean of 66 days. Snowmaking days in the Shanty Creek Basin were only slightly lower, with a mean of 91 in cool winters and 61 in warm. Despite the lower amount of suitable snowmaking days based on temperature, more snowmaking is done at Shanty Creek Resorts in warm winters compared to cool winters (Figure 5.4). Mean total snowmaking hours in cool winters at Shanty Creek was 747 using over 195,000 m3 of water compared to a mean total of 710 hours in warm winters and over 158,000 m3 of water (Table 5.1). While the increased snowmaking appears counterintuitive when considering the increased natural snow depths in cool winters, there are likely indirect influences leading to more snowmaking such as more open runs and higher visitor amounts during these snowier, cooler winters (Table 5.1). This is evidenced by the much earlier end to snowmaking in these cooler winters which typically stops around February compared to snowmaking in warm winters which persists throughout the season and into the spring months of March and April. 119 Figure 5.3. Winter type depth and snowmaking suitability. Panel a) is a violin plot of the distribution of days with the necessary snow depth for skiing in all ski basins in warm and cool winters. The wider the plot the more data points fall within that range. Horizontal lines denote the 25th, 50th and 75th quantiles. Warm winters generally have less of these days. Panel b) shows the number of days in ski basins in warm and cool winters with adequate temperatures for snowmaking. In warm winters there are less of these days. 120 Figure 5.4. Snowmaking water use and hours. a) Cumulative daily water use for snowmaking at Shanty Creek Resorts for each warm and cool winter. b) Cumulative snowmaking hours at Shanty Creek in those same winters. The vertical lines represent the holiday period. Cooler winters end with higher snowmaking hours and water used, but snowmaking ends several months earlier in the season than warm winters. As winter temperatures continue to warm in the future under current climate change projections, the adaptation strategies by ski slopes will likely focus on snowmaking. However, the increased pressures on natural resources such as energy and water necessary for snowmaking 121 will make this strategy less and less viable in the long term (Knowles, 2019). Reduced natural snow depths in the early winter season, combined with the inflexible dates of the economically critical holiday season, will increase pressure for larger snowmaking capabilities, requiring larger volumes of water (Wilson et al., 2018). Under high emissions scenarios, snowmaking needs in the Great Lakes are likely to double by midcentury and triple by the end of the century (Scott et al., 2021). This may be offset somewhat by the shift to higher surface water flows in the winter months instead of the spring due to earlier melting, but it’s unlikely to offset the increased snowmaking needs from reduced natural snow (Ford et al., 2020; 2021). Ultimately this is likely to lead to profit margins that are too thin for many ski areas to remain financially viable in the long-term future, causing the significant reduction of the amount of slopes in operation concluded by Scott et al. (2021). Increased winter temperatures and the subsequent reduction in snowfall will negatively impact the ski industry as the number of visitors and thus ticket sales declines, potentially leading to billions of dollars of lost revenue annually (Chapagain et al., 2018; Shih et al., 2009; Wobus et al., 2017). Some ski resorts are already planning for this eventuality and trying to adapt by utilizing more efficient snowmaking technologies and increasing non-winter recreational opportunities such as mountain biking (Knowles, 2019; Ooi, 2017). However, it does not appear that many in the ski industry are taking the problem seriously enough with their long-term planning as they haven’t begun to feel a negative enough effects and there remains many uncertainties about the future (Tashman and Rivera, 2015). There are a variety of influences driving this slow response, but the level of commitment to climate adaptation and mitigation strategies utilized by the industry as a whole is dependent on the level of risk they feel (Rivera and Clement, 2019). Too little risk and the problem is deemed unworthy of the 122 expenditures, and too much risk leads to the belief that adaptation measures are too costly for the lack of return, leaving a “goldilocks” zone of assessed risk that’s able to pressure ski firms and management to action. The economic risk from these winter climate changes is not just to the resorts themselves, but to the regions around the slopes. Ski tourism leads to tourism spending in other areas such as lodging (only 28% of skiers in Michigan stayed at the resort lodging with the rest staying at other hotels in the region or personal residences), travel costs (40% of Michigan ski visitors came from out of state in 2000), and food (46% of ski visits in Michigan were overnight visits) (Shih et al., 2009). Changes to the ski industry also effect the housing markets of the surrounding region, with higher average winter temperatures decreasing the housing prices in the area as home buyers will no longer be able to enjoy the winter recreational opportunities (Galinato and Tantihkarnchana, 2018). All of these potential economic losses, both direct and indirect, have led to increased risk assessment during financial negotiations with ski areas by real estate managers, banks and investors (Knowles, 2019). While the general economic forecast for the industry is poor, there are several assumptions in this study that contribute to uncertainty. First, the ski slope data utilized was limited. Detailed data from Shanty Creek covers only 7 years, leading to an inability to determine long term trends, and this is data from just one resort out of over 40 in the state. Actual data regarding costs and visitor amounts are generalized and regionally aggregated, limited its ability to be used for full economic quantification. This is in addition to the plethora of influential spending factors not related to snow amounts on the slopes such as the amount of disposable income available for the average visitor due to the overall health of the economy. 123 4. Conclusions Warming winter temperatures resulting from climate change are leading to decreased snowpack thicknesses throughout the winter, shortened snow season length and increased number of bare ground days during the winter. These changes to the seasonal snowpack are likely to lead to severe economic losses by the ski industry in Michigan and the rest of the country. While technological improvements to snowmaking efficiency have been able to offset the decreased snow depth thus far, they are unlikely to overcome the increased pressures as winter temperatures continue to rise. This is likely to lead to significant financial losses to the industry and the closure of many ski areas. These economic losses will not be limited to just the industry but are likely to be felt throughout the surrounding regions that rely on the tourism dollars the ski industry generates. Without enhanced adaptation and mitigation techniques in the industry and significant reduction in climate changing emissions, there is a challenging long- term future of the ski industry in Michigan. Acknowledgments This study would not have been possible without the data graciously provided by Shanty Creek Resorts CEO Pete Bigford. This research was supported by USDA-NIFA grant 2015- 68007-2313, NOAA grant NA12OAR4320071, and NASA grant NNX11AC72G along with Michigan State University’s Environmental Science and Public Policy Program and the Department of Earth and Environmental Sciences. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of USDA NIFA, NASA, or NOAA. 124 REFERENCES 125 REFERENCES Adam, J.C., Hamlet, A.F., Lettenmaier, D.P., 2009. Implications of global climate change for snowmelt hydrology in the twenty-first century, in: Hydrological Processes. pp. 962–972. https://doi.org/10.1002/hyp.7201 Anderson, J.R., Hardy, E.E., Roach, J.T., and Witmer, R.E., 1976. A land use and land cover classification system for use with remote sensor data. Geological Survey Professional Paper 964. Andresen, J.A., 2012. Historical Climate Trends in Michigan and the Great Lakes Region, in: Dietz, T., Bidwell, D. (Eds.), Climate Change in the Great Lakes Region: Navigating an Uncertain Future. Michigan State University Press, East Lansing, Michigan, pp. 17–34. Angel, J.R., Kunkel, K.E., 2010. The response of Great Lakes water levels to future climate scenarios with an emphasis on Lake Michigan-Huron. J. Great Lakes Res. 36, 51–58. https://doi.org/10.1016/j.jglr.2009.09.006 Argyilan, E.P., Forman, S.L., 2003. Lake Level Response to Seasonal Climatic Variability in the Lake Michigan-Huron System from 1920 to 1995. J. Great Lakes Res. 29, 488–500. https://doi.org/10.1016/S0380-1330(03)70453-5 Arora, V.K., Boer, G.J., 2001. Effects of simulated climate change on the hydrology of major river basins. J. Geophys. Res. 106, 3335–3348. https://doi.org/10.1029/2000JD900620 Avanzi, F., De Michele, C., Morin, S., Carmagnola, C.M., Ghezzi, A., Lejeune, Y., 2016. Model complexity and data requirements in snow hydrology: seeking a balance in practical applications. Hydrol. Process. 30, 2106–2118. https://doi.org/10.1002/hyp.10782 Barnett, T.P., Adam, J.C., Lettenmaier, D.P., 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 438, 303–309. https://doi.org/10.1038/nature04141 Barrett, A.P., 2003. National Operational Hydrologic Remote Sensing Center Snow Data Assimilation System (SNODAS) Products at NSIDC. NSIDC Spec. Rep. 11 19. Berghuijs, W.R., Woods, R. a, Hrachowitz, M., 2014. A precipitation shift from snow towards rain leads to a decrease in streamflow-supplement. Nat. Clim. Chang. 4, 583–586. https://doi.org/10.1038/NCLIMATE2246 Blöschl, G., Hall, J., Viglione, A., Perdigão, R.A.P., Parajka, J., Merz, B., Lun, D., Arheimer, B., Aronica, G.T., Bilibashi, A., Boháč, M., Bonacci, O., Borga, M., Čanjevac, I., Castellarin, A., Chirico, G.B., Claps, P., Frolova, N., Ganora, D., Gorbachova, L., Gül, A., Hannaford, 126 J., Harrigan, S., Kireeva, M., Kiss, A., Kjeldsen, T.R., Kohnová, S., Koskela, J.J., Ledvinka, O., Macdonald, N., Mavrova-Guirguinova, M., Mediero, L., Merz, R., Molnar, P., Montanari, A., Murphy, C., Osuch, M., Ovcharuk, V., Radevski, I., Salinas, J.L., Sauquet, E., Šraj, M., Szolgay, J., Volpi, E., Wilson, D., Zaimi, K., Živković, N., 2019. Changing climate both increases and decreases European river floods. Nature. https://doi.org/10.1038/s41586-019-1495-6 Boyer, C., Chaumont, D., Chartier, I., Roy, A.G., 2010. Impact of climate change on the hydrology of St. Lawrence tributaries. J. Hydrol. 384, 65–83. https://doi.org/10.1016/j.jhydrol.2010.01.011 Brubaker, K.L., Rango, A., 1996. Response of snowmelt hydrology to climate change. Water. Air. Soil Pollut. https://doi.org/10.1007/BF00619293 Burakowski, E.A., Wake, C.P., Braswell, B., Brown, D.P., 2008. Trends in wintertime climate in the northeastern United States: 1965-2005. J. Geophys. Res. Atmos. 113. https://doi.org/10.1029/2008JD009870 Burakowski, E., Magnusson, M., 2012. Climate Impacts on the Winter Tourism Economy in the United States. Prep. Prot. Our Winters Nat. Resour. Def. Counc. https://www.nrdc.org/sites/default/files/climate-i. https://doi.org/10.1038/ngeo2304 Burnett, A.W., Kirby, M.E., Mullins, H.T., Patterson, W.P., 2003. Increasing Great Lake-effect snowfall during the twentieth century: A regional response to global warming? J. Clim. 16, 3535–3542. https://doi.org/10.1175/1520-0442(2003)016<3535:IGLSDT>2.0.CO;2 Byun, K., Hamlet, A.F., 2018. Projected changes in future climate over the Midwest and Great Lakes region using downscaled CMIP5 ensembles. Int. J. Climatol. 38, e531–e553. https://doi.org/10.1002/joc.5388 Campbell, J.L., Driscoll, C.T., Pourmokhtarian, A., Hayhoe, K., 2011. Streamflow responses to past and projected future changes in climate at the Hubbard Brook Experimental Forest, New Hampshire, United States. Water Resour. Res. 47. https://doi.org/10.1029/2010WR009438 Chamberlain, S., Anderson, B., Salmon, M., Erickson, A., Potter, N., Stachelek, J., Simmons, A., Ram, K., Edmund, H., rOpenSci, 2020. rnoaa: ‘NOAA’ Weather Data from R. R Package version 0.8.5.9421. https://CRAN.R-project.org/package=rnoaa Champagne, O., Arain, M.A., Coulibaly, P., 2019. Atmospheric circulation amplifies shift of winter streamflow in southern Ontario. J. Hydrol. 578. https://doi.org/10.1016/j.jhydrol.2019.124051 Chapagain, B.P., Poudyal, N.C., Bowker, J.M., Askew, A.E., English, D.B.K., Hodges, D.G., 2018. Potential Effects of Climate on Downhill Skiing and Snowboarding Demand and 127 Value at U.S. National Forests. J. Park Recreat. Admi. 36, 75–96. https://doi.org/10.18666/jpra-2018-v36-i2-8365 Cherkauer, K.A., Sinha, T., 2010. Hydrologic impacts of projected future climate change in the Lake Michigan region. J. Great Lakes Res. 36, 33–50. https://doi.org/10.1016/j.jglr.2009.11.012 Chin, N., Byun, K., Hamlet, A.F., Cherkauer, K.A., 2018. Assessing potential winter weather response to climate change and implications for tourism in the U.S. Great Lakes and Midwest. J. Hydrol. Reg. Stud. 19, 42–56. https://doi.org/10.1016/j.ejrh.2018.06.005 Chin, N., Byun, K., Hamlet, A.F., Cherkauer, K.A., 2018. Assessing potential winter weather response to climate change and implications for tourism in the U.S. Great Lakes and Midwest. J. Hydrol. Reg. Stud. 19, 42–56. https://doi.org/10.1016/j.ejrh.2018.06.005 Cline, D.W., 1997. Effect of Seasonality of Snow Accumulation and Melt on Snow Surface Energy Exchanges at a Continental Alpine Site. J. Appl. Meteorol. 36, 32–51. https://doi.org/10.1175/1520-0450(1997)036<0032:EOSOSA>2.0.CO;2 Clow, D.W., 2010. Changes in the timing of snowmelt and streamflow in Colorado: A response to recent warming. J. Clim. 23, 2293–2306. https://doi.org/10.1175/2009JCLI2951.1 Clow, D.W., Nanus, L., Verdin, K.L., Schmidt, J., 2012. Evaluation of SNODAS snow depth and snow water equivalent estimates for the Colorado Rocky Mountains, USA. Hydrol. Process. 26, 2583–2591. https://doi.org/10.1002/hyp.9385 Conover, W.J., 1971. Practical Nonparametric Statistics. New York: John Wiley & Sons, 309- 314. Costa, D., Zhang, H., Levison, J., 2021. Impacts of climate change on groundwater in the Great Lakes Basin: A review. J. Great Lakes Res. 47, 1613–1625. https://doi.org/10.1016/j.jglr.2021.10.011 DeCicco, L., Hirsch, R., Lorenz, D., Read, J., Walker, J., Carr, L., Watkins, D., 2020. dataRetrieval: Retrieval Functions for USGS and EPA Hydrologic Water Quality Data. R Package version 2.7.5. https://CRAN.R-project.org/package=dataRetrieval Demaria, E.M.C., Palmer, R.N., Roundy, J.K., 2016. Regional climate change projections of streamflow characteristics in the Northeast and Midwest U.S. J. Hydrol. Reg. Stud. 5, 309– 323. https://doi.org/10.1016/j.ejrh.2015.11.007 Demiroglu, O.C., Turp, M.T., Ozturk, T., Kurnaz, M.L., 2016. Impact of climate change on natural snow reliability, snowmaking capacities, and wind conditions of ski resorts in northeast Turkey: A dynamical downscaling approach. Atmosphere (Basel). 7. https://doi.org/10.3390/atmos7040052 128 Dyer, J.L., Mote, T.L., 2006. Spatial variability and trends in observed snow depth over North America. Geophys. Res. Lett. 33. https://doi.org/10.1029/2006GL027258 Farrand, W., Bell, D., 1982. Quaternary Geology of Michigan. Feng, S., Hu, Q., 2007. Changes in winter snowfall/precipitation ratio in the contiguous United States. J. Geophys. Res. Atmos. 112. https://doi.org/10.1029/2007JD008397 Feinstein, D.T., Hunt, R.J., and Reeves, H.W., 2010. Regional groundwater-flow model of the Lake Michigan Basin in support of Great Lakes Basin water availability and use studies. Scientific Investigations Report 2010-5109. https://doi.org/10.3133/sir20105109 Ford, C.M., Kendall, A.D., Hyndman, D.W., 2020. Effects of shifting snowmelt regimes on the hydrology of non-alpine temperate landscapes. J. Hydrol. 590, 125517. https://doi.org/10.1016/j.jhydrol.2020.125517 Ford, C.M., Kendall, A.D., Hyndman, D.W., 2021. Snowpacks decrease and streamflows shift across the eastern US as winters warm. Sci. Total Environ. 793, 148483. https://doi.org/10.1016/j.scitotenv.2021.148483 Gan, T.Y., Barry, R.G., Gizaw, M., Gobena, A., Balaji, R., 2013. Changes in North American snowpacks for 1979-2007 detected from the snow water equivalent data of SMMR and SSM/I passive microwave and related climatic factors. J. Geophys. Res. Atmos. 118, 7682– 7697. https://doi.org/10.1002/jgrd.50507 Hagenstad, M., Burakowski, E., Hill, R., 2018. The Economic Contributions of Winter Sports in a Changing Climate. Prot. Our Winter POW 1–69. Hamilton, L.C., Lemcke-Stampone, M., Grimm, C., 2018. Cold winters warming? perceptions of climate change in the north country. Weather. Clim. Soc. 10, 641–652. https://doi.org/10.1175/WCAS-D-18-0020.1 Hamlet, A.F., Mote, P.W., Clark, M.P., Lettenmaier, D.P., 2005. Effects of temperature and precipitation variability on snowpack trends in the western United States. J. Clim. 18, 4545– 4561. https://doi.org/10.1175/JCLI3538.1 Hayhoe, K., VanDorn, J., Croley, T., Schlegal, N., Wuebbles, D., 2010. Regional climate change projections for Chicago and the US Great Lakes. J. Great Lakes Res. 36, 7–21. https://doi.org/10.1016/j.jglr.2010.03.012 Hayhoe, K., Wake, C.P., Huntington, T.G., Luo, L., Schwartz, M.D., Sheffield, J., Wood, E., Anderson, B., Bradbury, J., DeGaetano, A., Troy, T.J., Wolfe, D., 2007. Past and future changes in climate and hydrological indicators in the US Northeast. Clim. Dyn. 28, 381– 407. https://doi.org/10.1007/s00382-006-0187-8 129 Hedrick, A., Marshall, H.P., Winstral, A., Elder, K., Yueh, S., Cline, D., 2015. Independent evaluation of the SNODAS snow depth product using regional-scale lidar-derived measurements. Cryosphere 9, 13–23. https://doi.org/10.5194/tc-9-13-2015 Hidalgo, H.G., Das, T., Dettinger, M.D., Cayan, D.R., Pierce, D.W., Barnett, T.P., Bala, G., Mirin, A., Wood, A.W., Bonfils, C., Santer, B.D., Nozawa, T., 2009. Detection and Attribution of Streamflow Timing Changes to Climate Change in the Western United States. J. Clim. 22, 3838–3855. https://doi.org/10.1175/2009JCLI2470.1 Hodgkins, G.A., Dudley, R.W., Huntington, T.G., 2003. Changes in the timing of high river flows in New England over the 20th Century. J. Hydrol. 278, 244–252. https://doi.org/10.1016/S0022-1694(03)00155-0 Hodgkins, G.A., Dudley, R.W., 2006a. Changes in late-winter snowpack depth, water equivalent, and density in Maine, 1926-2004, in: Hydrological Processes. pp. 741–751. https://doi.org/10.1002/hyp.6111 Hodgkins, G.A., Dudley, R.W., 2006b. Changes in the timing of winter-spring streamflows in eastern North America, 1913-2002. Geophys. Res. Lett. 33. https://doi.org/10.1029/2005GL025593 Hodgkins, G.A., Dudley, R.W., Aichele, S.S., 2007. Historical Changes in Precipitation and Streamflow in the U . S . Great Lakes Basin , 1915 – 2004 Scientific Investigations Report 2007 – 5118, Program. Hollander, M., Wolfe, D.A., 1973. Nonparametric Statistical Methods. New York: John Wiley & Sons, 68-75. Homer, C., Huang, C., Yang, L., Wylie, B., Coan, M., 2004. Development of a 2001 national land-cover database for the United States. Photogramm. Eng. Remote Sens. https://doi.org/10.14358/PERS.70.7.829 Huntington, T.G., Hodgkins, G.A., Keim, B.D., Dudley, R.W., 2004. Changes in the proportion of precipitation occurring as snow in New England (1949-2000). J. Clim. 17, 2626–2636. https://doi.org/10.1175/1520-0442(2004)017<2626:CITPOP>2.0.CO;2 Hyndman, D.W., Kendall, A.D., Welty, N.R.H., 2007. Evaluating temporal and spatial variations in recharge and streamflow using the integrated landscape hydrology model (ILHM), in: Geophysical Monograph Series. https://doi.org/10.1029/171GM11 Isard, S.A., Schaetzl, R.J., 1998. Effects of winter weather conditions on soil freezing in southern Michigan. Phys. Geogr. 19, 71–94. Ishida, K., Ohara, N., Ercan, A., Jang, S., Trinh, T., Kavvas, M.L., Carr, K., Anderson, M.L., 2019. Impacts of climate change on snow accumulation and melting processes over 130 mountainous regions in Northern California during the 21st century. Sci. Total Environ. 685, 104–115. https://doi.org/10.1016/j.scitotenv.2019.05.255 Iwata, Y., Nemoto, M., Hasegawa, S., Yanai, Y., Kuwao, K., Hirota, T., 2011. Influence of rain, air temperature, and snow cover on subsequent spring-snowmelt infiltration into thin frozen soil layer in northern Japan. J. Hydrol. 401, 165–176. https://doi.org/10.1016/j.jhydrol.2011.02.019 Javed, A., Cheng, V.Y.S., Arhonditsis, G.B., 2019. Detection of spatial and temporal hydro- meteorological trends in Lake Michigan, Lake Huron and Georgian Bay. Aquat. Ecosyst. Heal. Manag. 22, 1–14. https://doi.org/10.1080/14634988.2018.1500850 Jefferson, A., Nolin, A., Lewis, S., Tague, C., 2008. Hydrogeologic controls on streamflow sensitivity to climate variation. Hydrol. Process. 22, 4371–4385. https://doi.org/10.1002/hyp.7041 Jennings, K.S., Winchell, T.S., Livneh, B., Molotch, N.P., 2018. Spatial variation of the rain- snow temperature threshold across the Northern Hemisphere. Nat. Commun. 9, 1–9. https://doi.org/10.1038/s41467-018-03629-7 Johnson, S.L., Stefan, H.G., 2006. Indicators of climate warming in Minnesota: Lake ice covers and snowmelt runoff. Clim. Change 75, 421–453. https://doi.org/10.1007/s10584-006-0356- 0 Johnston, C.A., Shmagin, B.A., 2008. Regionalization, seasonality, and trends of streamflow in the US Great Lakes Basin. J. Hydrol. 362, 69–88. https://doi.org/10.1016/j.jhydrol.2008.08.010 Kanno, Y., Letcher, B.H., Hitt, N.P., Boughton, D.A., Wofford, J.E.B., Zipkin, E.F., 2015. Seasonal weather patterns drive population vital rates and persistence in a stream fish. Glob. Chang. Biol. https://doi.org/10.1111/gcb.12837 Kendall, A.D., 2009. Predicting the Impacts of Land Use and Climate on Regional-Scale Hydrologic Fluxes. Michigan State University. L. B. Knowles, N., 2019. Can the North American ski industry attain climate resiliency? A modified Delphi survey on transformations towards sustainable tourism. J. Sustain. Tour. 27, 380–397. https://doi.org/10.1080/09669582.2019.1585440 Li, D., Wrzesien, M.L., Durand, M., Adam, J., Lettenmaier, D.P., 2017. How much runoff originates as snow in the western United States, and how will that change in the future? Geophys. Res. Lett. 44, 6163–6172. https://doi.org/10.1002/2017GL073551 Lilliefors, H.W., 1967. On the Kolmogorov-Smirnov Test for Normality with Mean and Variance Unknown. J. Am. Stat. Assoc. https://doi.org/10.1080/01621459.1967.10482916 131 Mankin, J.S., Diffenbaugh, N.S., 2015. Influence of temperature and precipitation variability on near-term snow trends. Clim. Dyn. 45, 1099–1116. https://doi.org/10.1007/s00382-014- 2357-4 McCabe, G.J., Wolock, D.M., 2010. Long-term variability in Northern Hemisphere snow cover and associations with warmer winters. Clim. Change. https://doi.org/10.1007/s10584-009- 9675-2 McCabe, G.J., Wolock, D.M., 2010. Long-term variability in Northern Hemisphere snow cover and associations with warmer winters. Clim. Change 99, 141–153. https://doi.org/10.1007/s10584-009-9675-2 MDEQ, 2008. WelLogic, statewide ground water database Menne, M.J., Durre, I., Korzeniewski, B., McNeal, S., Thomas, K., Yin, X., Anthony, S., Ray, R., Vose, R.S., E.Gleason, B., Houston, T.G., 2012. Global Historical Climatology Network - Daily (GHCN-Daily), Version 3 [WWW Document]. NOAA Natl. Clim. Data Cent. https://doi.org/10.7289/V5D21VHZ Michigan Department of Environment, Great Lakes and Energy (EGLE), 2019. Michigan Quaternary Geology. Available at http://gis- michigan.opendata.arcgis.com/datasets/quaternary-geology-features Mortsch, L., Hengeveld, H., Lister, M., Wenger, L., Lofgren, B., Quinn, F., Slivitzky, M., 2000. Climate Change Impacts on the Hydrology of the Great Lakes-St. Lawrence System. Can. Water Resour. J. 25, 153–179. https://doi.org/10.4296/cwrj2502153 Mote, P.W., 2003. Trends in snow water equivalent in the Pacific Northwest and their climatic causes. Geophys. Res. Lett. https://doi.org/10.1029/2003GL017258 Motiee, H., Mcbean, E., 2017. Assessment of Climate Change Impacts on Groundwater Recharge for Different Soil Types-Guelph Region in Grand River Basin, Canada. Ecopersia 5, 1731–1744. https://doi.org/10.18869/modares.Ecopersia.5.2.1731 Musselman, K.N., Addor, N., Vano, J.A., Molotch, N.P., 2021. Winter melt trends portend widespread declines in snow water resources. Nat. Clim. Chang. 11, 418–424. https://doi.org/10.1038/s41558-021-01014-9 Musselman, K.N., Clark, M.P., Liu, C., Ikeda, K., Rasmussen, R., 2017. Slower snowmelt in a warmer world. Nat. Clim. Chang. https://doi.org/10.1038/nclimate3225 Myneni, R., Knyazikhin, Y., and Park, T., 2015. MOD15A2H MODIS Leaf Area Index/FPAR 8- Day L4 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAC. http://doi.org/10.5067/MODIS/MOD15A2H.006 132 National Operational Hydrologic Remote Sensing Center (NOHRSC), 2004. Snow Data Assimilation System (SNODAS) Data Products at NSIDC, Version 1. [Snow Water Equivalent, 1 October 2003 to 31 May 2017]. Boulder, Colorado USA. NSIDC: National Snow and Ice Data Center. Doi: Https://doi.org/10.7265/N5TB14TC. [25 September 2017]. Naz, B.S., Kao, S.C., Ashfaq, M., Rastogi, D., Mei, R., Bowling, L.C., 2016. Regional hydrologic response to climate change in the conterminous United States using high- resolution hydroclimate simulations. Glob. Planet. Change 143, 100–117. https://doi.org/10.1016/j.gloplacha.2016.06.003 Notaro, M., Lorenz, D., Hoving, C., Schummer, M., 2014. Twenty-first-century projections of snowfall and winter severity across central-eastern North America. J. Clim. https://doi.org/10.1175/JCLI-D-13-00520.1 Novotny, E. V., Stefan, H.G., 2007. Stream flow in Minnesota: Indicator of climate change. J. Hydrol. 334, 319–333. https://doi.org/10.1016/j.jhydrol.2006.10.011 Nygren, M., Giese, M., Barthel, R., 2021. Recent trends in hydroclimate and groundwater levels in a region with seasonal frost cover. J. Hydrol. 602, 126732. https://doi.org/10.1016/j.jhydrol.2021.126732 Ooi, N., 2012. 8 Managing for resilience in the face of climate change: The adaptive capacity of U.S. ski areas. Peacock, S., 2012. Projected twenty-first-century changes in temperature, precipitation, and snow cover over north america in CCSM4. J. Clim. 25, 4405–4429. https://doi.org/10.1175/JCLI-D-11-00214.1 Perry, E., Manning, R., Xiao, X., Valliere, W., Reigner, N., 2018. Social Climate Change: The Advancing Extirpation of Snowmobilers in Vermont. J. Park Recreat. Admi. 36, 31–51. https://doi.org/10.18666/jpra-2018-v36-i2-8307 Pradhanang, S.M., Anandhi, A., Mukundan, R., Zion, M.S., Pierson, D.C., Schneiderman, E.M., Matonse, A., Frei, A., 2011. Application of SWAT model to assess snowpack development and streamflow in the Cannonsville watershed, New York, USA. Hydrol. Process. 25, 3268–3277. https://doi.org/10.1002/hyp.8171 Prism Climate Group, Oregon State University. 2018. Retrieved July 9, 2018, from http://prism.oregonstate.edu R Development Core Team, 2019. R: A Language and Environment for Statistical Computing. R Found. Stat. Comput. https://doi.org/10.1007/978-3-540-74686-7 Rivera, J., Clement, V., 2019. Business adaptation to climate change: American ski resorts and warmer temperatures. Bus. Strateg. Environ. 28, 1285–1301. https://doi.org/10.1002/bse.2316 133 Royston, P., 1995. Remarik AS R94: A remark on Algorithm AS 181: The W test for normality. Applied Statistics 44, 547-551. Scott, D., Steiger, R., Knowles, N., Fang, Y., 2020. Regional ski tourism risk to climate change: An inter-comparison of Eastern Canada and US Northeast markets. J. Sustain. Tour. 28, 568–586. https://doi.org/10.1080/09669582.2019.1684932 Scott, D., Steiger, R., Rutty, M., Knowles, N., Rushton, B., 2021. Future climate change risk in the US Midwestern ski industry. Tour. Manag. Perspect. 40, 100875. https://doi.org/10.1016/j.tmp.2021.100875 Seaber, P.R., Kapinos, F.P., Knapp, G.L., 1987. Hydrologic Unit Maps (USA). US Geol. Surv. Water-Supply Pap. 2294. https://doi.org/10.3133/wsp2294 Shanty Creek Resort, 2021. 2021 winter media kit winter wonderland. Shih, C., Holecek, D.F., 2020. Ticket Sales. Math. Teach. Middle Sch. 10, 130–139. https://doi.org/10.5951/mtms.10.3.0130 Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture. Soil Survey Geographic (SSURGO) Database for CONUS. Available online. Accessed January 25, 2020. Stewart, I.T., Cayan, D.R., Dettinger, M.D., 2004. Changes in snowmelt runoff timing in western North America under a `business as usual’ climate change scenario. Clim. Chang. 62, 217– 232. https://doi.org/10.1023/B:CLIM.0000013702.22656.e8 Stottlemyer, R., Toczydlowski, D., 2006. Effect of reduced winter precipitation and increased temperature on watershed solute flux, 1988-2002, Northern Michigan. Biogeochemistry 77, 409–440. https://doi.org/10.1007/s10533-005-1810-1 Sturm, M., Goldstein, M.A., Parr, C., 2017. Water Resources Research editor. Water Resour. Res. 53, 3534–3544. https://doi.org/10.1002/2017WR020840 Suriano, Z.J., Robinson, D.A., Leathers, D.J., 2019. Changing snow depth in the Great Lakes basin (USA): Implications and trends. Anthropocene 26, 100208. https://doi.org/10.1016/j.ancene.2019.100208 Tarboton, D.G., Luce, C.H., 1996. Utah energy balance snow accumulation and melt model (UEB). Tashman, P., Rivera, J., 2016. Ecological Uncertainty, Adaptation, and Mitigation in the U.S. Ski Resort Industry: Managing Resource Dependence and Institutional Pressures. Strateg. Manag. J. 37, 1507–1525. https://doi.org/10.1002/smj 134 Tu, J., 2009. Combined impact of climate and land use changes on streamflow and water quality in eastern Massachusetts, USA. J. Hydrol. 379, 268–283. https://doi.org/10.1016/j.jhydrol.2009.10.009 United States Geological Survey (USGS), 2011. NLCD 2011 Land Cover Coterminous United States. Available at https://mrlc.gov/data United States Geological Survey (USGS), 2014. Watershed Boundary Dataset (ver. USGS Watershed Boundary Dataset: Subbasin (8-digit) 4th level for the entire United States (published 20140311). United States Geological Survey (USGS), 2018. National Water Information System data available on the World Wide Web (USGS Water Data for the Nation). (n.d.). http://doi.org/http://dx.doi.org/10.5066/F7P55KJN United States Geological Survey (USGS), 2019. National Hydrography Dataset (NHD) available on the World Wide Web (USGS National Hydrography Products). https://www.usgs.gov/national-hydrography/national-hydrography-dataset United States Geological Survey (USGS), 2020. National Elevations Dataset available on the World Wide Web (USGS National Map Data). https://www.usgs.gov/the-national-map- data-delivery/gis-data-download U.S. Fish & Wildlife Service (2018). National Wetlands Inventory. U.S. Fish & Wildlife Service. https://data.nal.usda.gov/dataset/national-wetlands-inventory. Wilson, G., Green, M., Mack, K., 2018. Historical climate warming in the white mountains of New Hampshire (USA): Implications for snowmaking water needs at Ski Areas. Mt. Res. Dev. 38, 164–171. https://doi.org/10.1659/MRD-JOURNAL-D-17-00117 Wobus, C., Small, E.E., Hosterman, H., Mills, D., Stein, J., Rissing, M., Jones, R., Duckworth, M., Hall, R., Kolian, M., Creason, J., Martinich, J., 2017. Projected climate change impacts on skiing and snowmobiling: A case study of the United States. Glob. Environ. Chang. 45, 1–14. https://doi.org/10.1016/j.gloenvcha.2017.04.006 Wolock, D.M., Winter, T.C., McMahon, G., 2004. Delineation and Evaluation of Hydrologic- Landscape Regions in the United States Using Geographic Information System Tools and Multivariate Statistical Analyses. Environ. Manage. https://doi.org/10.1007/s00267-003- 5077-9 Xia, Y., Mitchell, K., Ek, M., Sheffield, J., Cosgrove, B., Wood, E., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., Lettenmaier, D., Koren, V., Duan, Q., Mo, K., Fan, Y., Mocko, D., 2012. Continental-scale water and energy flux analysis and validation for the North American Land Data Assimilation System project phase 2 (NLDAS-2): 1. Intercomparison and application of model products. J. Geophys. Res., 117, D03109, doi:10.1029/2011JD016048 135 Zou, B., Rockne, K.J., Vitousek, S., Noruzoliaee, M., 2018. Ecosystem and transportation infrastructure resilience in the great lakes. Environment 60, 18–31. https://doi.org/10.1080/00139157.2018.1495508 136