Northern forest winters have lost cold, snowy conditions that are important for ecosystems and human communities

Abstract Winter is an understudied but key period for the socioecological systems of northeastern North American forests. A growing awareness of the importance of the winter season to forest ecosystems and surrounding communities has inspired several decades of research, both across the northern forest and at other mid‐ and high‐latitude ecosystems around the globe. Despite these efforts, we lack a synthetic understanding of how winter climate change may impact hydrological and biogeochemical processes and the social and economic activities they support. Here, we take advantage of 100 years of meteorological observations across the northern forest region of the northeastern United States and eastern Canada to develop a suite of indicators that enable a cross‐cutting understanding of (1) how winter temperatures and snow cover have been changing and (2) how these shifts may impact both ecosystems and surrounding human communities. We show that cold and snow covered conditions have generally decreased over the past 100 years. These trends suggest positive outcomes for tree health as related to reduced fine root mortality and nutrient loss associated with winter frost but negative outcomes as related to the northward advancement and proliferation of forest insect pests. In addition to effects on vegetation, reductions in cold temperatures and snow cover are likely to have negative impacts on the ecology of the northern forest through impacts on water, soils, and wildlife. The overall loss of coldness and snow cover may also have negative consequences for logging and forest products, vector‐borne diseases, and human health, recreation, and tourism, and cultural practices, which together represent important social and economic dimensions for the northern forest region. These findings advance our understanding of how our changing winters may transform the socioecological system of a region that has been defined by the contrasting rhythm of the seasons. Our research also identifies a trajectory of change that informs our expectations for the future as the climate continues to warm.

missing more than five days of data within a single winter (e.g., Feng and Hu, 2007), ten days across a winter (e.g., Huntington et al., 2004;Feng and Hu, 2007;Kunkel et al., 2009), and fewer than 50% of stations over a ten-year time period (e.g., Huntington et al., 2004;Feng and Hu, 2007;Vincent et al., 2015). However, using such stringent criteria would have omitted nearly all of the Canadian sites from the analysis, thereby limiting the geographic scope of the study.
Because the period of record completeness differed among measurements recorded at any given site, a station could pass completeness criteria for one variable and fail for another. For example, total daily precipitation might pass, while snowfall might fail. Due to the dearth of stations in the USHCN and the NCDAEC databases that fit within our study area, we did not completely remove a site from the analysis unless both the precipitation and the temperature variables tested within a station failed the completeness criteria. Overall, 37 stations were retained for development and testing of winter climate change indicators. All 37 stations passed record completeness criteria for both precipitation and temperature, such that all 37 stations were included when evaluating winter climate change indicators. Average record completeness for each station was as follows: 79% for precipitation; 74% for daily maximum temperature, 69% for daily minimum temperature, 62% for snowfall, and 27% for snow depth, where record completeness was calculated as the total number of years containing winter records divided by the total number of years spanned by the record for each individual site by variable combination (Table S1). These values were somewhat lower than those reported in prior studies that used stricter completeness criteria (e.g., Huntington et al., 2004). Missing records were distributed evenly across the 100-year timeframe from 1916 to 2017 and were not more frequent earlier in the time series.

Section S2. Modeling Snow Depth and Snowfall
Snow depth and snowfall data were estimated for the entire dataset using a degree-day snowmelt model (Kokkonen et al., 2006;Buttle, 2009;Crossman et al., 2016) implemented in R (R Core Team, 2017) using the snow.sim function within the hydromad package (Andrews et al., 2011).
This model requires two driving variables: daily precipitation (mm d -1 ) and air temperature (°C).
The model simulates snowfall by assuming that above a threshold temperature (Tmax), all precipitation falls as rain, and below a threshold temperature (Tmin) it falls as snow, where Tmax must be greater than Tmin. It also simulates a snowmelt rate (kd in mm d -1 ) as occurring above a threshold temperature (where Tmelt = Tmin). Temperature thresholds and snowmelt rates were allowed to vary among different stations. Additional parameters can include a degree day factor for snow freezing (kf), liquid storage capacity of the snowpack (rcap), and correction factors for rain (cr) and snow (cs) to account for precipitation gage undercatch. These latter parameters (kf, rcap, cr, cs) were set to zero. The model can perform equally well without them in simulating snow water equivalent (SWE), and their omission reduces the overall complexity of the calibration step (Kokkonen et al., 2006). Implementation of the model required calibration of three parameters: a rainfall temperature threshold (Tmax), a snowfall temperature threshold (Tmin), and a melt rate (kd). When selecting parameter values, Tmax and Tmin were allowed to vary between -2 and 2°C, and kd ranged from 1 to 4 mm d -1 based on prior values published in the literature (e.g., Kokkonen et al., 2006). Since data from weather stations were reported as integers, the step increment for all parameter values was 1. Because Tmax has to be greater than Tmin for the model to function, this resulted in 50 possible unique parameter combinations of Tmax , Tmin, and kd for calibrating snow.sim. Prior implementation of the model has contained greater numbers of parameter combinations (e.g., Kokkonen et al., 2006). However, such simulations have included step increments as small as 0.05 for data collected to at least two significant digits, which would not be appropriate for our data set. For each parameter suite, the accumulation of the snowpack as SWE was simulated using the daily mass balance of snowfall minus snowmelt. The modeled snowpack that best matched the data for each weather station was selected by correlating the 50 possible simulated snow depths (in SWE) against the measured snow depth values (in solid snow; e.g., Crossman et al., 2016) and selecting the modeled SWE with the highest r 2 value, which varied from a minimum of 0.05 to a maximum of 0.80 with an average of 0.50. Snow depth values that had been recorded as zero in the historical record were removed prior to the correlation analysis to remove potential bias associated with observer error. Snowfall and rainfall were then simulated using the parameter set of Tmax and Tmin that had resulted in the best modeled SWE, as snow.sim only outputs total precipitation regardless of whether it falls as rain or snow.
Because measured and modeled snow depth and snowfall were not 1:1 analogs (measured were in solid snow, modeled were in SWE), calculations of winter climate change indicators from modeled snow data consisted of presence-absence metrics rather than shifts in amounts of snow cover or snowfall. The model accurately simulated the presence of a snowpack 98% of the time and precipitation falling as snow 95% of the time.