Projected changes in cold hardiness zones and suitable overwinter ranges of perennial crops over the United States

Average annual absolute minimum temperatures (TNn) provide a means of delineating agriculturally relevant climate zones and are used to define cold hardiness zones (CHZ) by the United States Department of Agriculture. Projected changes in TNn, mean winter minimum temperatures, and CHZs over the conterminous United States (CONUS) were assessed using an ensemble of statistically downscaled daily climate model output through the mid 21st century (2041–2070). Warming of TNn is on average ∼40% greater than that of mean winter minimum temperatures across CONUS with an average climate velocity of 21.4 km decade−1 resulting in widespread shifts in CHZs. These changes enable a geographic expansion of thermally suitable areas for the cultivation of cold-intolerant perennial agriculture including almond, kiwi, and orange crops. Beyond these crops, warming of TNn has broad implications for food security and biotic interactions.


Introduction
A number of studies have highlighted the importance of both mean and extreme minimum temperatures to ecological systems. For example, monthly average minimum temperatures have been used for habitat mapping (e.g. Ledig et al 2010), crop yield assessment (e.g. Lobell and Field 2007), and pest monitoring (e.g. Trần et al 2007, Paradis et al 2008. While the distribution of species may not be directly linked to mean annual or monthly temperatures, extreme minimum temperatures have established links to the overwinter survival rates of insects (e.g. Bale 1996, Stahl et al 2006 and plants (Alden andHermann 1971, Vetaas 2002). Consequently, extreme minimum temperatures provide constraints on the potential geographic range of natural and cultivated species (e.g. Woodward et al 2004), and can impact crop yields (e.g. Porter andGawith 1999, Gu et al 2008).
Cold damage to plants may occur at a range of minimum temperatures (T min ) depending on species sensitivity and phenological stage (Sakai andLarcher 1987, Larcher 2005). Given the challenges in generalizing plant cold tolerance, the United States Department of Agriculture (USDA) developed cold hardiness zones (CHZ) based on annual minimum T min (TN n ) averaged over a climatological period. These zones provide guidance for where plants might survive through winter temperatures and are presently used to establish nursery crop insurance standards. A map of CHZs was first published by the USDA in 1960, and was updated in 1990(Cathey 1990  Improved understanding of projected changes in temperature extremes-including TN n -have implications for informing climate adaptation approaches for crop cultivation, identifying areas at risk for invasive species expansions, and tracking potential changes in electricity and infrastructure needs. Increases in TN n under climate change will result in significant redistribution of biologically relevant thermoclines and subsequent species (e.g. Diffenbaugh et al 2008). Observed warming in annual average T min across the conterminous United States (CONUS) has resulted in poleward and altitudinal shifts in thermoclines, with spatially varying climate velocity-that is the speed and direction that a given property migrates with climate change (e.g. We build on these aforementioned studies by examining TN n and CHZ projections using an ensemble of global climate model output downscaled to a spatial resolution congruent with contemporary agroclimatic information, evaluating projected changes in TN n relative to mean winter (December-February) T min (TN , DJF ) and calculating the climate velocity of CHZs, TN n and TN . DJF Further, we complement previous work by examining the impact of projected changes in TN n on thermally suitable areas for the cultivation of three high market-value perennial fruit and nut crops: Nonpareil almond, Hayward kiwi, and Navel orange.

Data and methods
We obtained daily T min data from twenty global climate models (GCMs) participating in the fifth phase of the Climate Model Intercomparison Project (CMIP5) (Taylor et al 2012) that were statistically downscaled over CONUS using Multivariate Adaptive Constructed Analogs (MACA) method (Abatzoglou and Brown 2012) for both historical  and future  experiments. Downscaled data were trained using the gridded surface meteorological dataset of Abatzoglou (2013) at a 1/24th degree resolution grid that ensures that quantiles of the downscaled historical GCM period adhere to those of the observed record . The gridded dataset of Abatzoglou (2013) is a hybrid product that bias correct data from the North American Land Data Assimilation System (NLDAS2; Mitchell et al 2004) with monthly data from the Parameter Regression on Independent Slopes Model (PRISM; Daly et al 1994), and exhibits nominal biases for temperature extremes such as TN n when compared to in situ weather stations (i.e., coldest 1% of daily TN had a mean bias of +0.5°C compared to data from Global Historical Climate Network stations). MACA uses an analog approach for mapping GCM fields to observed fields and applies an equidistant quantile mapping bias correction procedure (Li et al 2010, Pierce et al 2015) that preserves the differences between future and historical daily temperatures from GCM simulations across quantiles, including TN n and other extreme values.
Dynamical downscaling using regional climate models (RCM) is arguably better suited for assessing climate extremes modulated by mesoscale land-surface phenomena (e.g., snow-albedo feedback). However, the restricted availability of RCM output from multiple GCMs and the additional statistical bias correction procedures needed for local assessment limited our analysis to the statistically downscaled products. We conduct a complementary analysis to facilitate a comparison between statistically downscaled products used in our analysis and dynamically downscaled results from two RCMs (CanCM4 and RCM4) using a common GCM ensemble member from the second generation Canadian Earth System Model (CanESM2) forced with RCP 8.5 as part of the CORDEX project (Giorgi et al 2009).
We constrained our analysis to model simulations for the historical period  and mid 21st century period (2041-2070). We chose to assess midcentury projections in TN n , TN , DJF and CHZs because of the limited ability for developing meaningful management strategies relevant to end-of-century projections. We primarily focus on future experiments run under Representative Concentration Pathway 8.5 (RCP 8.5) given that inter-model variability exceeds inter-scenario variability for these time horizons (Kharin et al 2007(Kharin et al , 2013, and emissions trajectories to date have more closely followed RCP 8.5 (Peters et al 2013).
TN n for each winter-centric year was calculated from November-March, along with TN . DJF We calculated 30 year averages of TN n and TN DJF for each model for both the historic and mid 21st century time periods and considered both multi-model ensemble averages, as well as the ensemble 25th and 75th percentiles to assess intermodel variability. The climate velocities of multimodel mean TN n and TN DJF between historical and mid 21st century were calculated using a distance-based velocity algorithm (Hamann et al 2015). This algorithm determines the shortest distance between locations with analog climates and divides by the number of years between the two climate periods to provide the climate velocity in units of km yr −1 . We calculated both forward (i.e. current-to-future) and backward (future-to-current) velocities of the ensemble average temperature and report the minimum of the two velocities as a conservative estimate.
The multi-model mean of 30 year average TN n values were also used to define CHZs. Hardiness zones range from −51.1°C to 21.1°C with each zone spanning 5.6°C and being comprised of half-zones A and B, each covering a 2.8°C range. Projected changes in CHZs and the velocity of CHZ shifts were calculated. While CHZ projections may be useful for assessing climate change impacts on crop cultivation, we utilize minimum temperature thresholds (TN CROP ) for dormancy as a means of examining how projected changes in TN n may expand thermally suitable areas for crop survival. We chose to examine the impacts of projected changes in TN n on Nonpareil almonds, Hayward kiwis, and Navel oranges because of their relatively high market value. These cultivars also provide examples across a range of hardiness threshold temperatures, from −25°C for Nonpareil almonds (Janick and Moore 1996), to −12°C for Hayward kiwifruit (Strik 2005), to −4.4°C for Navel oranges (Fake and Norton 2012). Using the multi-model mean TN n for both the historical and mid-century periods, we calculated the percent area over CONUS with TN n values above TN CROP . Additionally, to provide a more conservative measure of potential changes in crop cultivation area, we assessed the percent suitable land area where at least 80% of the models showed TN n >TN CROP .

Results
Ensemble average projected increases in TN DJF range from 1.7°C in the southeastern US to more than 5°C in the Upper Midwest and northern Great Plains ( figure 1(a)). While the spatial pattern of warming for TN DJF resembled that seen in TN n , the magnitude of warming of the latter was more acute across a majority of CONUS ( figure 1(b)). The ensemble average TN n warming ranged from 1.8°C to more than 7°C warming, yielding a 40% greater increase compared to TN DJF when averaged over CONUS. This results in an additional 2°C of warming of TN n over TN DJF across a broad region of the Midwestern US, Great Lakes and interior northwestern US ( figure 1(c)). This asymmetric warming was found for all downscaled GCMs across much of the northern half of the United States from the Great Plains to the Atlantic Ocean, as well as for much of the Intermountain West. Conversely, fewer GCMs showed differential warming across portions of the southern United States, the Rocky Mountains, and portions of the Southwest including California and Arizona (figure S1). Intermodel variability, represented by 25th and 75th percentiles of projected increases in TN n and TN DJF across models (figure S2), was largest over the northern US and the northern Rocky Mountains of Idaho and Montana, for TN DJF and TN n , respectively.
The velocity of TN n also varied spatially ( figure 2(a)). The mean (median) estimate of the speed of TN n over CONUS was 21.4 km decade −1 (16.2 km decade −1 ), albeit with substantial spatial heteorogeneity as seen in prior assessments of climate velocity (e.g., Loarie et al 2009, Dobrowski et al 2013). The fastest speed of TN n was found over the northern Great Plains and Midwestern US due to large increases in TN n coupled with a weak spatial gradient in TN n , while slow speeds were found along the West Coast, in the Southwest, and in coastal Florida. By comparison, the velocity of projected TN DJF was less than TN n , with a mean (median) of 15.6 km decade −1 (12 km decade −1 ) and with similar spatial patterns ( figure 2(b)). As an artifact of the spatial bounds of our data, forward-looking climate velocities have no analog climates within CONUS for parts of the Northern Plains and more localized areas in the Rocky Mountains and the Northeast (figure S3). Backward-looking velocities show analog climates over 95% of locations and differ from forward-looking velocities, particularly over the topographically complex Western US.
As CHZs are calculated from average TN n , those locations showing the largest warming of TN n also exhibited the largest projected increases in CHZs (e.g. from zone 5 to 6). CHZs of downscaled GCM data from historical runs (figure 3(a)) were similar to published CHZs from observational records (e.g. Daly et al 2012). A comparison of ensemble mean TN n downscaled from historical runs and TN n calculated using daily PRISM data from 1981 to 2010 showed absolute biases <1.5°C over CONUS, with a mean bias of +0.1°C, suggesting reasonable agreement. Mid-century CHZ projections showed northward and upward shifts in existing zones (figure 3(b)), with a mean (median) shift over CONUS of ∼93 km (∼56 km) by the mid 21st century. Nearly all (98%) of CONUS exhibited an increase in CHZ (i.e. toward warmer absolute minimum temperature) using the multi-model mean, and no location saw a decrease in CHZ. Similar changes in CHZ were projected by the mid 21st century using RCP 4.5 forcing (figure S4).
Warming TN n (and consequent shifting CHZs) resulted in an increase in land area with sufficiently warm temperatures for overwinter survival of crops. Over the historical period, approximately 24% and 5% of CONUS had sufficiently warm TN n for overwinter survival of oranges and kiwifruit, respectively (figures 4(a) and (c)). Mid 21st century projections of TN n would enable an expansion of land with suitable overwinter temperatures to approximately 37% and 9% of CONUS for kiwifruit and oranges, respectively; the extent of TN CROP for oranges expanded northward along coastal areas and kiwi expanded northward from its historical range (figures 4(b) and (d)). The majority (∼74%) of CONUS showed multi-model mean TN n >TN CROP for almonds over the historical period (figure 4(e)), expanding into the north central plains and covering ∼93% of CONUS by the mid 21st century (figure 4(f)). A more conservative approach, where at least 80% of the models have TN n >TN CROP , shows comparable results: the percent land area suitable for crop survival over the historical (future) period was ∼73% (∼90%) for almonds, ∼23% (∼32%) for kiwi, and ∼5% (∼8%) for oranges.
Similar patterns of warming are evident across the statistically and dynamically downscaled data, however changes in TN n and TN DJF are more heterogeneous in the dynamically downscaled outputs ( figure S5). The spatial correlation of changes in TN n (TN DJF ) between the downscaled data and the RCMs were 0.80 (0.87) for RCM4, and 0.83 (0.88) for RCA4. The raw GCM output, statistically downscaled data, and both RCMs show amplified warming of the TN n versus TN DJF over the majority of CONUS. Whereas the RCMs highlight heterogeneous warming in the topographically complex western United States, the inter-RCM variability is quite large.

Discussion and conclusions
The mechanisms responsible for the amplified warming of TN n are likely a function of Arctic amplification and land-atmosphere interactions. The Arctic and interior Canada are primary air mass source regions for cold air outbreaks over CONUS that typically result in TN n . Observed amplification in warming rates over high latitude landmasses and the poles versus the midlatitudes has contributed to an increase in the temperature of cold air masses that have impacted CONUS over the second half of the 20th century (Walsh et al 2001, Hankes andWalsh 2011). Huybers et al (2014) showed a pattern of decreased variance of the coldest 5% of TN in DJF with warming TN DJF on an interannual basis in observations, which supports the amplified warming of TN n . Continued amplified warming rates of source regions for cold air outbreaks likely contribute to the larger warming rate of TN n versus TN . DJF While changes in atmospheric circulation with climate change have been hypothesized to increase the potential for cold air outbreaks (e.g. Francis and Vavrus 2012), decreases in temperature variance as a result of climate change would reduce the potential for cold air outbreaks (e.g. Schneider et al 2015, Screen 2014). Changes in snow cover and depth can also increase warming rates as the high albedo and thermal emissivity of snow cover helps promote exceptionally cold temperatures. Consequently, projected declines in snowfall (Lute Climate velocity may shape the distribution of ecological zones and resident species (Loarie et al 2009). As TN n has a direct link to species viability, we suggest that the climate velocity of such metrics is important for changes in range shifts in agricultural and natural ecosystems. While their methodology for calculating climate velocities differs from that used here, Dobrowski et al (2013) showed similar patterns in the velocity of mean T min over the 20th century, though our mean projected velocities of TN n are greater than the average velocity of mean T min in that study. The greater velocity of TN n versus TN DJF suggests a hastened rate of change that may be important for planning and adaptation efforts, and in fact the velocity of change may be more important for some adaptation efforts than the magnitude of the change itself.
It is important to note the uncertainty in the projection of extremes in GCM data and associated statistical downscaling that may not fully capture the mesoscale land-surface feedbacks that can modify warming of temperature extremes. Statistically and dynamically downscaled products for a common GCM ensemble member generally show similar patterns in TN n and TN , DJF but localized magnitudes differ. The heterogeneity exhibited in the dynamically downscaled products (e.g. over the western US) likely highlight regions where snow-albedo feedbacks are captured by RCMs (e.g. Salathé et al 2008, Pepin et al 2015. In this respect, RCMs may help to unveil changes occurring at local scales that are not adequately resolved using statistical approaches. However, the influence of snow-albedo feedbacks is contingent upon accurately simulating snowcover changes. The differences in magnitude and spatial heterogeneity of the change in TN n and TN DJF between the RCMs examined here indicate the challenges in refining the magnitude of change at local scales. Likewise, the lack of ensemble GCM-RCM combinations and the potential for GCM biases to propagate into RCM simulations currently limit a comprehensive analysis suitable for research of this sort. However, coordinated experiments such as CORDEX (Giorgi et al 2009) can better elucidate uncertainty that arises through downscaling approaches as well as highlight value-added downscaling from RCMs on changes in TN n that may help refine our results.
The increase in TN n and subsequent shifts in CHZs projected for mid-century periods supports previous work on changing thermal suitability envelopes. For example, the analysis of Lobell et al (2006) on perennial crops over California showed favorability for future crop development at higher latitudes or elevations. Similarly, Olesen et al (2007) project an expansion in thermal suitability zones for maize production over Europe during the 21st century. While our analysis does not consider other factors (e.g. heat tolerance thresholds, chilling hour requirements, water availability, competing land use) that govern where crops can be cultivated, warming of TN n may provide opportunities for crop production in regions that are currently thermally limited by cold extremes. However, there are many caveats to the potential for crop expansion with respect to warming TN n . For the perennials examined here that are either early blooming or highly sensitive to frost damage, commercial cultivation occurs almost exclusively in areas where TN n is much warmer than TN CROP and there are few studies providing thorough examination of threshold temperatures for cold hardiness (Janick and Moore 1996). It should be noted that the TN CROP values used in this study are temperatures that would severely damage or kill crops during overwinter dormancy; during other phenological stages, crops may be at higher risk for damage from less extreme cold temperatures. Further, while these threshold temperatures may be tolerated during dormancy for a few hours, many hours below TN CROP would result in increased damage or mortality (Fake and Norton 2012). Additionally, tolerance may decrease on nights with little wind when radiative heat loss can cool plant tissues below the ambient air temperature (Johnson 2011).
Warming TN n and projected shifts in CHZs have implications for agricultural and natural vegetation, land management, the energy sector and infrastructure. In addition to cultivated crops, native and invasive species and pests may also see geographic expansion, resulting in additional challenges for agricultural land managers as well as those managing forests, rangelands and other natural resources (e.g. Noss 2001). Moreover, an increase in TN n may also have economic impacts. Provided that the greatest electrical demands for heating occur during the coldest temperatures, the anticipated reductions in heating demand assessed from projected changes in TN DJF may be augmented further with greater rates of warming of TN n (Scott andHuang 2007, Mideksa andKallbekken 2010). In addition to lowered heating costs, further economic impacts of warming TN n include the reduced cost of transportation infrastructure repairs as warmer T min extremes reduce thermal stress on asphalt and damage from frost heaves (e.g. Mills and Andrey 2002).
The differential warming exhibited between changes in mean and extreme minimum winter temperatures highlights the importance of assessing both means and extremes in understanding potential impacts of climate change. Through utilizing daily projections to illustrate results with direct implications for climate change impacts, we show the benefit in revisiting previous studies whose analyses were limited temporally and spatially by previously unavailable downscaled daily data, and suggest that for applied purposes statistically downscaled products may be preferable to RCMs for multi-member ensemble studies. Finally, although the caveats presented above highlight the need for additional research to more fully account for the role of climatological factors governing crop survival, our results show promise for geographic expansion of thermally limited cultivars under climate change.