Spring plant phenology and false springs in the conterminous US during the 21st century

The onset of spring plant growth has shifted earlier in the year over the past several decades due to rising global temperatures. Earlier spring onset may cause phenological mismatches between the availability of plant resources and dependent animals, and potentially lead to more false springs, when subsequent freezing temperatures damage new plant growth. We used the extended spring indices to project changes in spring onset, defined by leaf out and by first bloom, and predicted false springs until 2100 in the conterminous United States (US) using statistically-downscaled climate projections from the Coupled Model Intercomparison Project 5 ensemble. Averaged over our study region, the median shift in spring onset was 23 days earlier in the Representative Concentration Pathway 8.5 scenario with particularly large shifts in the Western US and the Great Plains. Spatial variation in phenology was due to the influence of short-term temperature changes around the time of spring onset versus season-long accumulation of warm temperatures. False spring risk increased in the Great Plains and portions of the Midwest, but remained constant or decreased elsewhere. We conclude that global climate change may have complex and spatially variable effects on spring onset and false springs, making local predictions of change difficult.


Introduction
The onset of spring plant growth, or 'spring onset,' has shifted earlier in the year in recent decades amid rising global temperatures (Cleland et al 2007, Ault et al 2011, McCabe et al 2012. Whereas a longer growing season may increase carbon uptake and potentially mitigate climate change (Black et al 2000, Dragoni et al 2011, earlier leaf and flower emergence has led to phenological mismatches between plant resources and many of those animals depending upon them (Walther et al 2002, Schweiger et al 2008, Saino et al 2011, Burkle et al 2013, Kellermann and van Riper IIIrd 2015. Earlier spring onset can also lead to increased risk of false springs, when subsequent hard freezes damage new, vulnerable plant growth in ecological and agricultural systems (Cannell and Smith 1986, Inouye 2008, Knudson 2012. At many locations, projected future increases in global temperatures will likely result in earlier spring onset and fewer frost days overall. However, the rates and magnitude of these changes, determining the likelihood of phenological mismatches and false springs, are not known.
Studies of spring plant phenology often define spring onset with one of two easily identified events: leaf emergence ('leaf out') or flower emergence ('first bloom'; Schwartz 1998, Wolfe et al 2005, Polgar and Primack 2011. The timing of these events are largely determined by temperature and photoperiod (Bernier 1988, Bernier et al 1993, Polgar and Primack 2011, but the exact phenological response varies among plant species, ecotypes, and genotypes (Schwartz 1993, Polgar andPrimack 2011). Despite this variation, general phenological models, like the widely-used spring indices (SI-x) (Schwartz 1990, Ault et al 2015, capture the behavior of a wide variety of plants in natural and agricultural systems (Wolfe et al 2005, Schwartz and Hanes 2010. Historical studies using the SI-x (Schwartz et al 2006, Ault et al 2011, McCabe et al 2012 and other general indices (Peterson and Abatzoglou 2014) indicate that spring onset has generally shifted earlier, but there is considerable variation in the magnitude of trends, and in some areas spring onset now occurs later than in the recent past. Both leaf out and first bloom will likely continue to shift earlier in the year with rising global temperatures, but the complex nature of this process makes it difficult to make a priori estimates of how changes in temperature will affect spring onset. However, newly available daily records from climate projections allow for high resolution modeling of spring onset that may help identify spatial patterns of change, even in topographically complex areas where temperatures may vary considerably over short distances.
Many plants are resistant to cold weather while dormant through the winter. However, sub-freezing temperatures after spring onset can damage vulnerable plant tissue, and reproductive growth stages later in spring typically make plants more susceptible to damage from cold (Sakai et al 1987, Augspurger 2013. Damage due to false springs is often observed in natural systems (Inouye 2000, 2008, Augspurger 2011, and lost plant productivity can negatively impact dependent animal populations (Blais et al 1955, Nixon andMcClain 1969). False springs can also strongly affect agricultural systems (Paulsen andHeyne 1983, Brown andBlackburn 1987). For example, the false spring of 2012 caused $500 million in damages to fruit and vegetables in Michigan (Knudson 2012. Broad-scale studies suggest that cold weather is diminishing more rapidly than changes in spring onset, therefore decreasing the risk of a false spring at most locations (Schwartz et al 2006, Marino et al 2011, Peterson and Abatzoglou 2014. However, several recent false springs in North America suggest that the risk remains (Gu et al 2008, Hufkens et al 2012, Knudson 2012 and has even increased in some areas (Inouye 2008, Augspurger 2013. Empirical evidence is mixed, and the future risk of false springs is still largely unknown.
In this study, we asked two questions: first, how will spring onset change by 2100 in the conterminous US, and second, how will relative changes in spring onset and freezing conditions affect the likelihood of future false springs?

Methods
Downscaled climate data We based our spring onset and false spring calculations on daily minimum and maximum temperature data for 1950-2100 from the Coupled Model Intercomparison Project 5 (CMIP5) multi-model ensemble General Circulation Models (GCM) dataset, statistically downscaled to 1/8th degree resolution with the bias-corrected constructed analog (BCCA) technique (Maurer et al 2007, Bureau of Reclamation 2014. With the BCCA technique, the GCM data are first debiased with historical records at the coarse resolution of the GCM and then downscaled to a 1/8th degree grid using linear combinations of past historical conditions (Maurer et al 2010). We downloaded data covering much of North America from the Global Organization for Earth System Science Portals (Maurer et al 2007) but restricted our analysis to the conterminous US. We analyzed BCCA-derived data from 19 GCMs and two emissions scenarios: Representative Concentration Pathways (RCPs) 4.5 (medium-low emissions) and 8.5 (high emissions) (table 1). For conciseness, we focus the main text on the model averages and on RCP8.5, though we do include some analysis of inter-model variability. We present results for RCP4.5 in appendix A, and extensive results for individual models under both scenarios in appendix B.
Extended SI-x We calculated spring onset using the extended SI-x metrics . The original spring index formulas were empirically derived to predict leaf out and first bloom of lilac (Syringa chinesis 'Red Rothomagensis') and two honeysuckle clones (Lonicera tatarica 'Arnold Red' and L. korolkowii 'Zabeli') in the US (Schwartz 1993). However, these formulas have proven useful as a general model of plant phenology, as the first leaf and first bloom dates are well correlated with those of many natural and agricultural plants (Wolfe et al 2005, Schwartz and Hanes 2010 and with the 'green up' of ecosystems worldwide (Schwartz 1990, Schwartz et al 2006. Additionally, first bloom captures leaf out for later spring species, like many trees (Schwartz et al 2006, Schwartz andHanes 2010). The SI-x models that we applied have been modified to estimate leaf out and first bloom in both temperate and subtropical environments , Ault et al 2015. We defined spring onset using both leaf out and first bloom The use of both definitions provided 'bookends' to capture the spring phenology of a wide variety of plant species.
The SI-x include an equation to predict leaf out for each of the three plant species using empiricallyderived equations based on the photoperiod (represented by the day of the year), short-term growing degree hours (GDHs), and the season-long cumulative count of high-energy synoptic events, hereafter 'warm spells.' GDHs are the number of hours above 0.6°C, summed over two three-day periods (current and 5-7 days earlier). Hourly temperatures were interpolated from daily minimum and maximum temperatures (Ault et al 2015). The season long warm spell count is defined as the cumulative number of 3 day periods with GDH>637. The formula for each species is updated for each day of year, and the first day of the year with a SI-x score exceeding a predetermined threshold is considered the date of leaf out (see Ault et al 2015 for full models and coefficients). The average day of the year in which this occurs among the lilac and two honeysuckle clones is considered the 'leaf out date' (Schwartz et al 2006. The first bloom calculations begin after leaf out for each species. First bloom depends on two factors: the number of days since leaf out and the accumulated GDHs (AGDH) since leaf out (Ault et al 2015). The base GDHs are calculated in the same manner as for leaf out, but the AGDH is cumulative rather than based only on recent temperatures. The date of first bloom occurs for each species when first bloom values reach 1000 (Ault et al 2015). Again, we used the average day of year of the first bloom of the three species as our first bloom date in subsequent analyses. Changes in first bloom can only result from changes in AGDH, the final value of which is relatively constant given the threshold score of 1000. Therefore, to measure changes in first bloom we examined the duration of time between leaf out and first bloom. Within each year, we considered dates through 31 July, which was sufficient to capture spring growth even at high elevations while avoiding the return of colder temperatures accompanying the onset of autumn.
We compared changes in spring onset from the CMIP5 simulated historical , mid-century , and end of the century (2071-2100) time periods. We summarized changes by Omernik ecoregions (Omernik 1987), i.e., areas with relatively similar topography and weather patterns, by comparing probability density functions for during the historical and end-century time periods. We also used these results to examine differences in inter-model variability among ecoregions. For each ecoregion, we calculated coefficient of variation among model density kernels at a given day of year, and then taking the average of these values. To prevent values near zeros from overriding the signal, we removed from consideration days of the year in which the mean density was less than 1/50th of the peak value for that ecoregion. To determine if changes were robust among models at each location in the presented maps, we took the mean of each spring onset date from each model during the historical period, and compared it to the same values from future time periods using a Welch's t-test to allow for unequal variances. Non-significant changes were masked with white in the map figures. Finally, we have included more extensive results from individual models in appendix B.

False springs
The exact temperature at which plant tissue is damaged depends upon the species and environmental conditions, but a daily minimum temperature below −2.2°C is likely to cause damage, and is considered a 'hard freeze' (Schwartz 1993, Marino et al 2011, Peterson and Abatzoglou 2014. Furthermore, the later the freeze occurs after plant growth begins, the more damage is likely to occur because plants are in a more susceptible phenological phase (Schwartz et al 2006, Marino et al 2011, Peterson and Abatzoglou 2014. Because we calculated two measures of spring plant phenology, we also had two measures of false springs. We defined an early false spring event as the occurrence of a hard freeze seven or more days after leaf out, following Peterson and Abatzoglou (2014). Similarly, we defined a late false spring event as a hard freeze any time after first bloom, when plants are more susceptible to freezing conditions. We summarized changes in false springs by ecoregion, and tested for changes that were significant among models, in the manner described for spring onset. To investigate spatial variation in false springs, we calculated daily minimum temperatures around the leaf out and first bloom dates at each location. For each year, we calculated the average daily minimum temperature for the 7 days before and after the date of interest, and averaged these annual values for each CMIP5 time period.

Data validation
We validated the leaf out, first bloom, and false springs values calculated from the BCCA downscaled climate models versus those calculated from Maurer gridded historical temperature records (Maurer et al 2002), because the latter were the basis for downscaling of the BCCA data itself (Bureau of Reclamation 2014). We compared mean values, interannual variability, and trends between Maurer and the model average BCCA values. The BCCA-derived values compared favorably with historical records (appendix C). To increase calculation speed, we wrote code in C# based on earlier Matlab documentation and code (Ault et al 2015), with a few small modifications to fit our study (see appendix D). The code is available from the spring index Github page.

Spring onset
During the historical period, model average spring onset, defined by leaf out, followed a latitudinal gradient in the Eastern half of the US, whereas elevation played a larger role in the West ( figure 1(a)).
Along the Gulf Coast and the Southern portion of the North American Desert region, spring onset was driven primarily by temperatures immediately preceding leaf out (i.e., GDH; figure 2(a)). Just North of these areas, season-long weather conditions were more important (i.e., cumulative warm spells), and the influence of warm spells was particularly strong along the West Coast and Pacific Northwest ( figure 2(b)). At the highest latitudes and elevations, temperatures were lower (figure 2) and photoperiod was most important (results not shown).
By the middle and end of the century, model average leaf out is projected to shift earlier at all locations and under both climate scenarios (RCP8.5 in figure 1; RCP4.5 in appendix A). For RCP8.5, the region-wide median change was 22.3 days by the end of the century. The largest changes in leaf out occurred in the Western US, in the North American Desert, Temperate Sierras, and Northwestern Forested Mountain ecoregions, with median shifts of 26.5-28.5 days earlier by the end of the century. The smallest shifts occurred in areas where leaf out was already early in the historical period (e.g., Gulf Coast, Desert Southwest; figure 1(e)). The largest changes in leaf out were caused by an increasing influence of cumulative warm spells (spatial correlation=0.87, figure 2(d)), particularly at high elevations and in the middle latitudes of the Great Plains and Eastern Temperate Forest ecoregions (figure 1(e)). Correspondingly, smaller shifts in leaf out occurred in areas where the importance of the GDH terms increased (spatial correlation=−0.66; figure 2(b)). An interesting exception was the Marine West Coast ecoregion, which saw a large advance in leaf out (figures 1(c) and (e)) while the influence of GDH increased in an area previously dominated by cumulative warm spells (figure 2).
The historical pattern of dates and projected changes in model average spring onset defined by first bloom were closely related to those of leaf out, partly because the calculation of first bloom begins at leaf out (figure 1). Region-wide changes in leaf out and first bloom were strongly correlated (ρ=0.86) and the RCP8.5 nationwide median change of 22.1 days by the end of the century for first bloom was very similar to that of leaf out, although the band of largest shifts in leaf out date is located North of the maximum shifts in first bloom (figures 1(e) and (f)). However, interesting spatial patterns emerge in the duration of the period between leaf out and first bloom (figure 3). The primary factor in triggering first bloom is the AGDH since leaf out (see methods). Therefore, in areas where change in leaf out was driven by GDH (e.g., the Deep South), AGDH increased rapidly as well, and thus the time between leaf out and first bloom decreased. This was particularly evident in the Marine West Coast ecoregion, and across the Southern Great Plains and Eastern Forest ecoregions ( figure 3(b)). However, we observed a different pattern in areas where leaf out was strongly accelerated by cumulative warm spell increases, such as the broad swath of maximum gain in warm spell points from the Northern Rockies Eastward to New England (figure 2(d)). Though both dates shifted earlier in the year (figures 1(e) and (f)), leaf out shifted more in these areas, and thus the time between leaf out and first bloom increased ( figure 3(b)). Though model average changes in both definitions of spring onset were significant between models at all locations (figures 1(c)-(f)), individual models did vary in their predictions. Individual model predictions of spring onset varied within years, but all models tracked the model mean trend (figure 4). When pooled by ecoregion, distributions of spring onset dates were shaped by geographical (e.g., latitudinal extent) and topographical features (figure 5). The shapes of these distributions changed by end-century, with earlier spring onset dates and more inter-model variability than seen in the historical time period (figure 5). However, levels of inter-model variability were still low (see coefficients of variation for each ecoregion in figure 5).

False springs
During the historical period, early false springs based on leaf out were very common throughout much of the region, occurring almost annually at many locations ( figure 6(a)). The only areas without frequent false springs based on leaf out were those that either rarely experienced freezing weather (e.g., the Gulf and West Coasts), or where leaf out was driven primarily by photoperiod (e.g., upper Midwest, Northern New England) and occurred later in the year ( figure 7(a)). Late false springs based on the first bloom occurred much less frequently, but they were most common in certain high elevation regions, as well as the Southern US ( figure 6(b)), where first bloom occurs early in the year ( figure 1(b)).
In future projections, false springs declined in frequency under both definitions and throughout much of the range (figures 6(c)-(f)) as temperatures generally moderated around the time of spring onset (figures 7(c) and (d)). The decline was especially prominent for early false springs because they were so     common during the historical period. However, large areas of the Great Plains and Eastern Forest ecoregions were projected to experience an increase in false springs, particularly for the late false springs that were relatively rare during the historical period (figures 6(d) and (f)). In these areas, leaf out and first bloom dates advanced so rapidly that simulated daily minimum temperatures at the time of spring onset were actually lower at the end of the century than during the historical period (figures 7(c) and (d), and demonstrated in figure 4(d)), thereby promoting the chances of a hard freeze.

Discussion
Our projections indicated that spring onset will occur earlier throughout the conterminous US by the year 2100. Spatial variation in the rates of change stemmed from differing contributions by projected temperature . We found similar variability in individual BCCA models ( figure 5, appendix B), but the year-to-year variability in these modelled atmospheric patterns is not expected to match historical records. Our averaging of many models canceled out the effects of the simulated internal variability and emphasized the role of greenhouse-forced climate change.
The extended SI-x that we calculated are well correlated with spring phenology in ecosystems worldwide , and the time from leaf out to first bloom spans the period of spring phenology for many plant species. However, the phenological response to environmental conditions varies among plant species and even among locations (Polgar and Primack 2011), and other phenological models have been developed. While the extended SI-x incorporate spring degree-above-threshold and photoperiod (day of year) measures to predict leaf out, winter chilling days are included in other models, reflecting the physiology of plant species that require cold temperatures to break dormancy (Polgar and Primack 2011). Our projections may not capture changes in spring phenology of such plants. For example, a study using a black ash (Fraxinus nigra) phenological model and a CMIP3 GCM predicted that spring leaf out will arrive later in the mid-latitudes of the Eastern US (Morin et al 2009), where we projected much earlier spring onset. Detailed phenological models may be required for species-specific models, whereas our study was designed to capture ecosystem green up in general, for which the extended SI-x are ideal .
Future changes in false spring risk depend on the relative change in the timing of spring onset and the last spring freeze (Cannell andSmith 1986, Peterson andAbatzoglou 2014). While both relate to changes in temperature, differences arise because spring phenology is a cumulative process throughout the spring whereas a single cold night can cause a hard freeze. Though our early definition was common enough to strain the term false spring, the shared spatial patterns indicate that future patterns of damaging false springs will likely be somewhere between our two 'bookend' scenarios ( figure 4). We projected only slightly earlier spring onset along the Gulf Coast, where projected increases in temperature nearly eliminated freezing temperatures, causing a large reduction in the risk of false springs. Similarly, decreases in false springs risk in the Northern forests were the result of a projected decrease in late cold temperatures, and only a slightly earlier spring onset in a region where photoperiod remained an important component of the phenology. This difference in timescales was particularly evident in the Great Plains ecoregion, where the cumulative warm spells caused large shifts in leaf out and first bloom. These dates changed more than the projected mean temperature increases would suggest, causing temperatures around those dates to become cooler at the end of the century than historically. Although it may at first glance be counterintuitive that increasing mean temperatures could increase the risk of a false spring, the phenomenon has already been reported in historical studies (Inouye 2008, Augspurger 2013. Plants face an evolutionary tradeoff between the benefit of earlier leaf emergence, and hence, a longer growing season, and the risk of tissue damage from a false spring. Variation in internal responses to environmental cues may allow some individuals to survive the opposing selective pressures from earlier leaf out and risk of tissue damage (Leinonen and Hanninen 2002, Gömöry and Paule 2011). Though extended SI-x are static, our results suggest these selective pressures will continue. Throughout most of the US a similar or reduced risk of false springs indicates that plants could continue their present responses to environmental cues, and the resulting earlier spring growth would have little negative consequence, but it may mean missing the benefit of a longer growing season. However in the Great Plains, conditions may favor individuals that demonstrate delayed spring leaf out and avoid the increased risk of tissue damage from false springs.
Given changing plant phenology, maintaining current plant-animal interactions will require evolutionary or behavioral adaptation by dependent animals (Visser andHolleman 2001, Berteaux et al 2004). However, animals may respond to different environmental cues than plants. For example, longdistance migratory birds respond to cues present in their overwintering habitat, such as day length, while plants in their summer breeding grounds respond to local environmental cues like temperature. Birds that have adapted to migrate earlier have maintained their population levels, while birds that retained historical temporal patterns in migrations have declined, at least in part due to phenological mismatches with plant-based resources (Saino et al 2011, Clausen andClausen 2013). Increasing temperatures have led to poor synchronization between moth emergence and leaf out of host trees (Visser and Holleman 2001). Ultimately, the ability of a species to respond to rapid phenological changes will depend upon generation time, levels of genetic variability, and the plasticity of phenological and behavioral traits (Berteaux et al 2004).
In summary, our projections indicate that widespread historical advances in spring plant phenology will continue into the future, albeit with considerable spatial variation in the rates of change and the risk of false springs. Extensive regional variation emphasizes the need for future predictions that are even more fine-scale and species specific, to better understand the potential effects on natural and agricultural systems. To facilitate such research, we have created an online repository of weather data. The data presented here and more can be downloaded at http://silvis.forest. wisc.edu.

Acknowledgments
We thank M Schwartz and T Ault for advice regarding the spring indices, S Schmidt for technical assistance, and B Bateman, J Gorzo, and R Behnke for helpful discussion of these ideas. Financial support was provided NASA Biodiversity Program and the Climate and Biological Response funding opportunity (NNH10ZDA001N-BIOCLIM). The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the US Fish and Wildlife Service. Any use of trade, product, or firm names are for descriptive purposes only and do not imply endorsement by the US Government. All associated data are available upon request or from http://silvis.forest.wisc.edu.

Appendix C. Data validation with Maurer historical data
Here, we conduct a validation of leaf out, first bloom, and false springs values calculated from the BCCA downscaled climate models versus those calculated from Maurer gridded historical temperature records (Maurer et al 2002). This dataset provided a natural comparison because it was used in the downscaling of the BCCA data itself (Bureau of Reclamation 2014).
The mean values of leaf out and first bloom were very similar between the Maurer and BCCA results (figures C1 and C2 top). Differences do exist, but they are small and not readily viewable on the map. Interannual variability in leaf out and first bloom were higher in the BCCA models, though we observed similar spatial patterns of variability with the Maurer data. There was considerable spatial variation in trends of spring onset through the base period for both Maurer (figures C1 and C2) and individual BCCA models (not shown). However, the BCCA model average smoothed over most of this variability, and there was typically a slight shift earlier in these dates through time. For both definitions of false springs, spatial patterns were nearly identical between the Maurer and BCCA model average values (figures C1, C2, bottom).  (table 1).

Appendix D. Comparison of C# code against Matlab package
To rapidly process our large set of data, we wrote code to calculate the extended SI-x and associated values in the C# language. We followed the documentation provided in Ault et al (2015) in writing the code. In this appendix we evaluate output from our code against the recent version of the associated Matlab script (si_ml_v6.0.0) associated with Ault et al (2015), and use test weather station data provided in the Matlab package. We found that in almost all instances, the C# code produced exactly the same results as the Matlab code for leaf out, first bloom, and all of the predictor variables therein. However, in rare cases, we found differences in leaf out date. First bloom was unaffected other than as a results of a change in leaf out date, so we do not consider it further here. After thoroughly analyzing the code, we found that the issue stems from a line in the Matlab code (Line 61 in leaf.m) that skips all processing on days that where the maximum temperature is below the baseline threshold. This was likely intended to save processing time on days that cannot contribute GDH, or to match legacy programs, but unfortunately this causes problems in lagged GDH values that the program stores (e.g., Line 124-5 in leaf.m).
We show an example below using the DDE2 component of SI-x, the sum of GDH over the current day and previous two days. In table D1, the GDH has not been calculated for the days of year 5 and 6 because the daily maximum temperature was below the baseline temperature. However, we can see that the lag variable has not been updated, and the DDE2 value for day of year 7 is the sum of days 3, 4, and 7 (highlighted in bold), rather than 5, 6, and 7. Values from the C# code or modified Matlab code also produce zero GDH for days 5 and 6, but these zeros are properly included in the DDE2 total for day 7 (table D2).
Because these values occur when temperatures are low, these changes in, e.g., DDE2, do not directly cause changes in the leaf out date. Instead, occasionally this type of error can cause the loss or gain of a warm spell, which are cumulative throughout the year (see methods). While this issue can shift the date of spring onset by a few days at a location within a year, the issue was rare and exploratory analysis showed no noticeable difference in the time period averages in the main text, within or among models, using the unmodified Matlab code or the C# code.