CLIMATOLOGY OF SPRINGTIME FREEZE EVENTS IN THE GREAT LAKES REGION AND THEIR IMPACT ON SOUR CHERRY YIELDS IN HISTORICAL AND PROJECTED FUTURE TIME FRAMES  By  Lydia Rill          A THESIS Submitted to  Michigan State University in partial fulfillment of the requirements for the degree of Geography - Master of Science 2016     ABSTRACT  CLIMATOLOGY OF SPRINGTIME FREEZE EVENTS IN THE GREAT LAKES REGION AND THEIR IMPACT ON SOUR CHERRY YIELDS IN HISTORICAL AND PROJECTED FUTURE TIME FRAMES  By  Lydia Rill Production of sour cherries has a significant impact on the economy of the Great Lakes Region, valued at more than $74 million. In contrast to cereal crops where water is the most limiting factor, perennials in temperate regions are limited by freeze injury, especially in the spring following initial phenological development, which was highlighted in 2002 and 2012 when yields decreased considerably. This study analyzed the spatial and temporal variability of springtime freeze events in the Great Lakes Region and their impact on sour cherry production, as well as explored the use of gridded climatic datasets. Additionally, this study examined the historical trends in sour cherry yield and potential future changes by the mid-century. The Great Lakes played a major role in the spatial variability of springtime freezes, as locations inland experienced colder temperatures than near the coast. Simulated damaging freeze events over the past 50 years were most common during the later phenological stages of the crop, while the most severe damage occurred in the earlier stages. Gridded datasets were less suitable for this application than individual station observations. Over time, phenological development has tended to begin earlier, and this trend was projected to continue into the mid-century at three stations located in northwest, west central, and southwest Lower Michigan. An ensemble of climate projections indicated a large uncertainty envelope surrounding changes in simulated sour cherry yield. The projected changes varied by emissions scenario, downscaling method, and climate model.  iii  ACKNOWLEDGEMENTS  I would like to thank my advisors, Jeff Andresen and Julie Winkler, and committee member, Roy Black, for their help and guidance in this research. They have been very supportive throughout graduate school, and I am very thankful for all that they have taught me. I would like to acknowledge the Department of Geography for its financial support. I enjoyed being part of the department; it expanded my interests and taught me a great deal. I would like to thank my family for all their support, specifically my parents, Brenda and Peter Rill, and my brother, Elliott. They are always understanding and supportive in everything I do.  I would also like to thank to my roommate, Deanna Apps, for all her help and support these last two years. She has been my go-to for my questions throughout graduate school and has kept me going through the late nights.    iv  TABLE OF CONTENTS  LIST OF TABLES ......................................................................................................................... vi LIST OF FIGURES ..................................................................................................................... viii KEY TO ABBREVIATIONS ........................................................................................................ xi CHAPTER 1: INTRODUCTION ................................................................................................... 1  CHAPTER 2: IMPACT OF SPRINGTIME FREEZE EVENTS ................................................... 7 2.1 Background ......................................................................................................................... 7 2.1.1 Vulnerability of Perennial Crops to Freezing Temperatures ..................................... 7 2.1.2 Springtime Freeze Types and their Impact on Perennial Crops ................................ 8 2.1.3 Influence of Phenology on Freeze Damage ............................................................. 10 2.1.4 Frost Risk of Perennial Crops .................................................................................. 11 2.1.5 Study Objectives ...................................................................................................... 14 2.2 Data and Methods ............................................................................................................. 15 2.2.1 Springtime Freeze Types ......................................................................................... 15 2.2.2 Sour Cherry Simulation Model ................................................................................ 17 2.2.3 Sour Cherry Simulations over Space ....................................................................... 18 2.3 Results and Discussion ..................................................................................................... 22 2.3.1 Temporal and Spatial Variability of Springtime Freezes ........................................ 22 2.3.2 Characteristics of Damaging Freezes....................................................................... 36 2.3.3 Influence of Various Climatic Input Data on Sour Cherry Simulations .................. 41  CHAPTER 3: HISTORICAL AND FUTURE YIELDS .............................................................. 49 3.1 Background ....................................................................................................................... 49 3.1.1 Susceptibility of Perennial Crops to Climate Change .............................................. 49 3.1.2 Historical Trends in Phenology and Frost Risk ....................................................... 50 3.1.3 Projected Future Changes in Phenology and Frost Risk .......................................... 52 3.1.4 Climate Projections .................................................................................................. 56 3.1.5 Climate Change Impacts on Perennial Crops in Michigan ...................................... 58 3.1.6 Study Objectives ...................................................................................................... 60 3.2 Data and Methods ............................................................................................................. 61 3.2.1 Sour Cherry Yield Model......................................................................................... 61 3.2.2 Historical Climate Data............................................................................................ 62 3.2.3 Sources of Future Climate Projections .................................................................... 62 3.2.4 Downscaling and Debiasing Methods ..................................................................... 66 3.2.5 Application of the Climate Projections .................................................................... 69 3.3 Results ............................................................................................................................... 70 3.3.1 Historical Temporal Trends ..................................................................................... 70 3.3.2 Projected Future Changes in Climatic Variables ..................................................... 77 3.3.3 Projected Future Changes in Sour Cherry Production ............................................. 88 3.3.3.1 Changes in Sour Cherry Yield ....................................................................... 88 v  3.3.3.2 Changes in Buds Remaining .......................................................................... 95 3.3.3.3 Changes in Poor Pollination Days .................................................................. 95 3.3.3.4 Changes in Phenology .................................................................................. 104 3.3.3.5 Changes in Frequency and Severity of Damaging Freeze Events ............... 113 3.4 Summary and Discussion ................................................................................................ 123  CHAPTER 4: CONCLUSION ................................................................................................... 128 APPENDIX ................................................................................................................................. 132 BIBLIOGRAPHY ....................................................................................................................... 134    vi  LIST OF TABLES  Table 2.1: Base 4ºC GDD Accumulations for Phenological Stages ............................................. 18  Table 2.2: Various Climatic Input Data for the Sour Cherry Simulation Model ......................... 20  Table 2.3: Frequency of Freezes at Traverse City, Muskegon, and South Bend, 1960-2015 ....... 23  Table 2.4: Frequency of Freezes with 2.5+ cm of Snow Cover at Traverse City, Muskegon, and South Bend, 1960-2015 ................................................................................................................. 24  Table 2.5: Frequency of Simulated Damaging Freezes at Traverse City, Muskegon, and South Bend, 1960-2015 ........................................................................................................................... 25  Table 2.6: Average Minimum Temperature during Radiation and Non-Radiation Freeze Events, 1960-2015 ..................................................................................................................................... 26  Table 2.7: Average Minimum Temperatures during Radiation Freeze Events with 2.5+ cm of Snow Cover, 1960-2015 ................................................................................................................ 31  Table 2.8: Average Minimum Temperatures during Simulated Damaging Freeze Events, 1960-2015............................................................................................................................................... 34  Table 2.9: Linear Regression for Observed Yields and Model Output ......................................... 42  Table 3.1: Global Climate Models used in the CMIP5 Projections ............................................. 64  Table 3.2: Climate Models used in the NARCCAP Projections ................................................... 65  Table 3.3: Available NARCCAP Simulations ............................................................................... 65  Table 3.4: Downscaling and Debiasing Methods ......................................................................... 67  Table 3.5: Linear Trends in Model Output at Eau Claire, Hart, and Maple City, 1960-2015..... 71  Table 3.6: Historical Averages of Yield Model Output, 1980-2000 ............................................. 76  Table 3.7: Multi-Model Mean Changes by the Mid-Century in Temperature and Precipitation . 78  Table 3.8: Simulated Changes (percent) by 2040-2060 in Sour Cherry Yield by Climate Projection Type ............................................................................................................................. 90  Table 3.9: Simulated Changes (percent) by 2040-2060 in Buds Remaining by Climate Projection Type ............................................................................................................................................... 96 vii   Table 3.10: Simulated Changes by 2040-2060 in Poor Pollination Days by Climate Projection Type ............................................................................................................................................. 103  Table 3.11: Simulated Changes by 2040-2060 in Date of Phenological Stage 2 (Side Green) by Climate Projection Type ............................................................................................................. 105  Table 3.12: Simulated Changes by 2040-2060 in Date of Phenological Stage 8 (Full Bloom) by Climate Projection Type ............................................................................................................. 109  Table 3.13: Simulated Changes by 2040-2060 in Frequency of Damaging Freezes by Climate Projection Type ........................................................................................................................... 115  Table 3.14: Simulated Changes (percent) by 2040-2060 in Average Severity of Damaging Freezes by Climate Projection Type ........................................................................................... 122   viii  LIST OF FIGURES  Figure 2.1: Study Weather Station Sites ........................................................................................ 16  Figure 2.2: Interpolated Average Minimum Temperature of March Freeze Events, 1960-2015 . 30  Figure 2.3: Interpolated Average Minimum Temperature of March Radiation Type Freeze Events, 1960-2015......................................................................................................................... 33  Figure 2.4: Hourly Characteristics of Simulated Damaging Radiation Type Freeze Events ....... 37  Figure 2.5: Hourly Characteristics of Simulated Damaging Non-radiation Type Freeze Events 38  Figure 2.6: Average Frequency and Severity of Simulated Damaging Freezes by Phenological Stage, 1960-2015 .......................................................................................................................... 40  Figure 2.7: Model Simulated Buds Remaining and Observed Yields for the Northwest Region .. 44  Figure 2.8: Model Simulated Buds Remaining and Observed Yields for the West Central Region....................................................................................................................................................... 46  Figure 3.1: Time Series of Simulated Yield, Buds Remaining, and Poor Pollination Days at Eau Claire, Hart, and Maple City, 1960-2015..................................................................................... 72  Figure 3.2: Time Series of Simulated Dates of Stage 2 and Stage 9 at Eau Claire, Hart, and Maple City, 1960-2015 ................................................................................................................. 74  Figure 3.3: Time Series of Simulated Damaging Freeze Events and Average Severity of Damage at Eau Claire, Hart, and Maple City, 1960-2015 ......................................................................... 75  Figure 3.4: Projected Change by 2040-2060 in Average Annual Maximum Temperature (°C) at Eau Claire by Climate Projection................................................................................................. 79  Figure 3.5: Projected Change by 2040-2060 in Average Annual Minimum Temperature (°C) at Eau Claire by Climate Projection................................................................................................. 80  Figure 3.6: Projected Change by 2040-2060 in Average Annual Precipitation per Wet Day (mm per day) at Eau Claire by Climate Projection .............................................................................. 82  Figure 3.7: Projected Change by 2040-2060 in Average Annual Frequency of Wet Days at Eau Claire by Climate Projection ........................................................................................................ 84  Figure 3.8: Monthly Deltas (°C) for Maximum Temperature from the CMIP5 Models ............... 85  ix  Figure 3.9: Monthly Deltas (°C) for Minimum Temperature from the CMIP5 Models ................ 86  Figure 3.10: Temperature Deltas for the NARCCAP Simulations ................................................ 87  Figure 3.11: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Eau Claire by Climate Projection ........................................................................................................ 89  Figure 3.12: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Hart by Climate Projection ........................................................................................................................ 92  Figure 3.13: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Maple City by Climate Projection............................................................................................................ 94  Figure 3.14: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Eau Claire by Climate Projection ........................................................................................................ 97  Figure 3.15: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Hart by Climate Projection ........................................................................................................................ 98  Figure 3.16: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Maple City by Climate Projection............................................................................................................ 99  Figure 3.17: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Eau Claire by Climate Projection............................................................................ 100  Figure 3.18: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Hart by Climate Projection ...................................................................................... 101  Figure 3.19: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Maple City by Climate Projection ........................................................................... 102  Figure 3.20: Simulated Change by 2040-2060 in Average Date of Stage 2 at Eau Claire by Climate Projection ...................................................................................................................... 106  Figure 3.21: Simulated Change by 2040-2060 in Average Date of Stage 2 at Hart by Climate Projection .................................................................................................................................... 107  Figure 3.22: Simulated Change by 2040-2060 in Average Date of Stage 2 at Maple City by Climate Projection ...................................................................................................................... 108  Figure 3.23: Simulated Change by 2040-2060 in Average Date of Stage 8 at Eau Claire by Climate Projection ...................................................................................................................... 110  Figure 3.24: Simulated Change by 2040-2060 in Average Date of Stage 8 at Hart by Climate Projection .................................................................................................................................... 111  x  Figure 3.25: Simulated Change by 2040-2060 in Average Date of Stage 8 at Maple City by Climate Projection ...................................................................................................................... 112  Figure 3.26: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Eau Claire by Climate Projection ........................................................... 114  Figure 3.27: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Hart by Climate Projection ..................................................................... 117  Figure 3.28: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Maple City by Climate Projection ........................................................... 118  Figure 3.29: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Eau Claire by Climate Projection ............................................................................. 119  Figure 3.30: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Hart by Climate Projection ....................................................................................... 120  Figure 3.31: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Maple City by Climate Projection ............................................................................. 121    xi  KEY TO ABBREVIATIONS  CMIP5  Coupled Model Intercomparison Project 5th phase COOP  Cooperative Observer Program GCM  Global Climate Model GDD  Growing Degree Day MLCR  Multiple Linear Contour Regression MLR  Multiple Linear Regression NARCCAP  North American Regional Climate Change Assessment Program NARR  North American Regional Reanalysis NASS  National Agricultural Statistics Service NLDAS-2  North American Land Data Assimilation System Phase 2 NCEI  National Centers for Environmental Information NCEP  National Center for Environmental Prediction NWS  National Weather Service PRISM  Parameter-elevation Regressions on Independent Slopes Model QM  Quantile Mapping RCM  Regional Climate Model RCP  Representative Concentration Pathway SRES  Special Report on Emissions Scenarios US  United States USDA  United States Department of Agriculture  1  CHAPTER 1: INTRODUCTION Perennial tree fruit production is an important segment of local and regional agricultural economies around the world. In the Great Lakes Region of the U.S.A., sour cherry is a major crop, with the majority of total national production in the State of Michigan (NASS, 2015). However, cherry yields vary greatly from year to year. For example, in 2002 and 2012 sour to springtime freeze events following unusually warm spells (NASS, 2002 and 2012). These dramatic losses inspired this research on sour cherry production in the Great Lakes Region, and the relationship between production and climate.  Perennial tree fruit yields are strongly influenced by the climate through springtime freezes, winter chill fulfillment, pollination, and plant disease (Winkler et al., 2013). In winter, cool temperatures are necessary for perennial trees to vernalize and blossom the following spring. However, in mild climates such as in California, the associated chilling requirement may not be sufficiently fulfilled, severely reducing yields (Luedeling et al., 2009). In colder, temperate climates, while the winter chill requirement is easily met, extreme cold temperatures during winter may limit over winter survival and the extent of the fruit growing regions. For example, in Canada winter injury due to extreme cold temperatures is a major constraint to apple production (Rochette et al., 2004). Cold temperatures can also cause damage in the fall before winter dormancy and in the spring after the trees break dormancy (Rodrigo, 2000).  Springtime freezes are in general the most limiting factor for fruit production in temperate climates (Rodrigo, 2000; Winkler et al., 2013). During the spring, freeze injury can occur in plants due to ice formation, specifically extra-cellular freezing (Rodrigo 2000; Snyder and De Melo-Abreu, 2005). As ice forms outside the cells, water within the cells evaporates and 2  deposits on the ice crystals, causing the ice to grow and the cells to desiccate, sometimes to the extent where the cells collapse and die (Snyder and De Melo-Abreu, 2005). Flowers and fruit can often recover and continue development after spring freezes, depending on the amount of damage to the vital tissues, but effects from the damage may be seen in subsequent fruit malformations (Rodrigo, 2000). To evaluate the risk of springtime freeze damage on perennials, phenology must be included, as fruit trees become increasingly vulnerable to freeze damage as they progress through various developmental stages (Chmielewski, 2013; Longstroth, 2005). Because environmental temperature is the main driver of phenology, warm spells in late winter or early spring can accelerate the rate of growth and development and cause the trees to be more susceptible to freeze damage (Vitasse et al., 2014). However, phenology is also dependent on the chilling requirement and photoperiod, so attempts to estimate or simulate phenology may become complex (Vitasse et al., 2014). Furthermore, the cultivar type of the perennial tree also affects the timing of growth (Chmielewski, 2013; Moghadam et al., 2009; San Martino et al., 2005). Although complex, accounting for phenology is crucial for quantifying the impact of springtime freeze events on fruit trees.  There are two main types of freezes, advection and radiation, that can damage perennial crops. Advection freezes occur when colder, sub-freezing air moves via wind into a region or area of interest, while radiation freezes occur under calm, clear conditions at night where the ground or plant canopy surface cools radiatively more quickly than the air above it. Microclimate directly influences the severity of the freezes; for example, topography and proximity to large water bodies can affect the temperatures (Winkler et al., 2013). In the Great Lakes Region, fruit tends to be grown within close proximity to the lakes, as they moderate the temperatures in areas 3  downwind of the lakes. Furthermore, orchards are typically placed in relatively high elevations because during radiation freeze events the cold air near the surface flows downslope and away from the trees (Perry, 1998). The type of freeze is critical, however, as the relatively high locations may cause the trees to be more exposed to wind during advection freezes and relatively more vulnerable to cold injury (Winkler et al., 2012b).   Scientific literature on perennial springtime freeze risk is limited and spread out geographically. Some studies use the date of the last frost with respect to the date of a phenological stage, usually bloom, as a proxy indicator of risk (Kurlus et al., 2013), while others include the frequency and severity of freezes in their research (Fitchett et al., 2014; Ladányi et al., 2010). Phenology observations are often used when available, but phenology may also be simulated, typically through the use of agro-climatic indices such as growing degree days (Winkler et al., 2013). In some studies, bud break is simulated using a chilling model (Molitor et al., 2014), and in others bud break is assumed to occur by a certain calendar date (Cittadini et al., 2006). Studies that quantify springtime freeze damage amounts are rare, as only one study in South America (Cittadini et al., 2006) and one in Turkey (Öztekin, 2008) estimated damage using critical temperature thresholds. There is a general lack of perennial crop simulation models as the slow growth of perennials tends to limit the amount of data for model development (Lobell and Field, 2011). Modelling cold injury risk for perennials is complex as the timing of bud break, phenology, and freezes are important. Additionally, studies on perennials should ideally account for the microclimate influences of the region of interest as environmental conditions may vary significantly over short distances. Chapter 2 of this study aims to better understand the spatial and physical characteristics of springtime freezes, and to quantify the effects of springtime 4  freezes on fruit production in the Great Lakes Region. Furthermore, Chapter 3 of this study examines the effects of a variable and changing climate on sour cherry yields in Michigan.  As the climate changes in the future, there is general uncertainty in how perennial crops will be impacted. Historically, there have been temporal trends in increasing temperatures and consequently earlier dates of phenology in many production areas (Chmielewski et al., 2004; Guédon and Legave, 2008; Kurlus et al., 2013; Legave et al., 2013). In the Midwest, the winters have been warming, allowing shifts in the spatial extents of plant hardiness regions (Daly et al., 2010). Additionally, the date of the last spring freeze has been advancing, but the interannual variability of cold outbreaks has been large in the past few decades (Kunkel et al., 2013). Although temperatures have been warming globally, historical trends in springtime freeze risk vary by location. For example, freeze risk increased in Illinois and Iran (Augspurger, 2013; Fitchett et al., 2014), and freeze risk decreased for perennial trees in eastern China and central Europe (Dai et al., 2013; Scheifinger et al., 2002). With temperatures projected to continue increasing in the future, phenology may advance earlier, resulting in an increased risk of frost damage. However, the timing and frequency of freezes may also shift, with little change or decreased risk (Fitchett et al., 2014; Winkler et al., 2013). Changes in climate variability may also significantly impact agriculture (Thornton et al., 2014), as temperature fluctuations can cause plants to be more susceptible to springtime freeze damage (Gu et al., 2008). Additionally, changes in temperature and precipitation in the future may influence pollination conditions for perennial crops. The rate of phenological advancement over time will play a large role in determining if springtime freeze risk will increase or decrease in the future. Many studies projected that phenology will advance to earlier dates in the future (Hur and Ahn, 2014). However, the results 5  for changes in future freeze risk vary. Some studies projected decreasing frost risk for fruit crops in Europe (Eccel et al., 2009; Hoffmann and Rath, 2013; Ladányi et al., 2010; Molitor et al., 2014), while other studies projected increasing frost risk (Chmielewski et al., 2010; Mosedale et al., 2015). Additionally, a couple of studies projected increasing frost risk in some regions and decreasing frost risk in others (Kaukoranta et al., 2010; Rochette et al., 2004). Furthermore, the future projections for freeze risk in Michigan exhibited high uncertainty (Winkler et al., 2012b). The inconsistent results may be due to differing climate projections as well as the methods of estimating freeze risk. Overall, the uncertainty surrounding frost risk for perennials in the future is large. Because the results of previous studies are inconsistent with each other, the use of multiple climate projections, creating an ensemble, helps to understand and quantify some of the uncertainty. Ensembles are important because various climate models, emissions scenarios, and downscaling methods can lead to different future projections (Hanssen-Bauer et al., 2005). Downscaling methods are commonly used as global climate models (GCMs) have coarse spatial resolutions and typically need to be downscaled to a region or location, especially for studies where regional or local influences have large impacts (Winkler et al., 2011a). Therefore for this study, a large ensemble of climate projections using varying emissions scenarios, climate models, and downscaling methods is used to evaluate possible future trends in sour cherry yields.  Overall a better understanding of the relationship between springtime freezes and fruit production in temperate climates is desired, especially for North America. Additionally, more information is needed on how fruit production may be impacted by global warming. In this study, Chapter 2 examines the historical relationship between springtime freezes and sour cherry phenology and damage in the Great Lakes Region for 1960-2015. The frequency and spatial 6  variability of damaging springtime freezes as well as the timing with respect to phenology are analyzed using a sour cherry model that simulates phenology and cold damage. Furthermore, the use of gridded datasets in place of climate observations in the sour cherry model is explored. In Chapter 3, the sour cherry model, with an additional yield algorithm accounting for springtime damaging freezes and pollination conditions, is used to determine historical and future trends in yield in the Great Lakes Region. A climate ensemble with 88 individual scenarios from various climate models, emissions scenarios, and downscaling methods is used to evaluate uncertainty in future changes in yield from a historical period, 1980-2000, to the mid-century, 2040-2060. 7  CHAPTER 2: IMPACT OF SPRINGTIME FREEZE EVENTS 2.1 Background 2.1.1 Vulnerability of Perennial Crops to Freezing Temperatures Sour cherries are a major component of the agricultural industry in Michigan and the Great Lakes Region (Winkler et al., 2013), with production in 2013 valued at more than $74 million in Michigan alone (MDARD, 2014). Michigan is the top producing state in the U.S. providing approximately 60-75% of the total sour cherry production during the past few years (NASS, 2015). The majority of production occurs in western Lower Michigan where regional climate is modified by the nearby Great Lakes and by topography. Production of temperate, perennial tree fruit crops in the region is significantly limited by freeze injury, in contrast to common cereal crops where water is the most limiting factor (Andresen et al., 2001; Charrier et al., 2015). Freezes during the spring season have a profound impact on sour cherry production in the region. For example, severe crop damage associated with spring freezes occurred in both 2002 and 2012. In 2012 production decreased 97% relative to the previous year, from 71 to 2 million kilograms due to record high temperatures during March of that year followed by a series of freezes in April and May (NASS, 2012). Similarly, sour cherry production in Michigan in 2002 dropped 95% from the previous year, from 135 to 7 million kilograms, due to unusually warm weather in April followed by a series of freeze events (NASS, 2002). Knowledge of the frequency and severity of freeze events is essential in quantifying the climatic-harm perennial crops in the fall before dormancy, during winter dormancy, and in the spring after the crop breaks dormancy and begins to grow (Rodrigo, 2000). During the winter, the tree is in endodormancy, hardened to withstand cold temperatures, but once chilling and photoperiod 8  requirements are fulfilled, the tree enters the ecodormancy stage, where warmer temperatures cause deacclimation, allowing the trees to become more susceptible to cold damage (Vitasse et al., 2014). While cold temperatures in winter can harm perennial trees, springtime freezes have a much larger effect on interannual variability of yield (Winkler et al., 2013). The timing of springtime freeze events based on the rate and stage of phenological development is a critical factor in freeze damage. Springtime freezes cause injury to the plant through extracellular ice formation, where the water inside the cells evaporates and deposits on the ice crystals, causing the ice to grow and the cells to desiccate, potentially leading to the collapse and death of cells (Snyder and De Melo-Abreu, 2005).  2.1.2 Springtime Freeze Types and their Impact on Perennial Crops The severity of springtime freeze damage to sour cherries and other perennial tree fruit crops depends on the magnitude of the temperature during the freeze, the duration of the freezing temperatures, and the timing of the freeze with respect to the phenological stage of the crop. The air temperatures during freezes in the Great Lakes Region are greatly influenced by microclimate; consequently, the locations of orchards tend to be in areas of modified climates, influenced by topography and proximity to large water bodies (Winkler et al., 2013). The two main categories of springtime freezes are radiation and advection freezes. Radiation freezes tend to occur under calm, clear weather conditions which allow temperatures to decrease radiatively during the overnight hours, typically resulting in the coldest temperatures at or near the ground surface and the formation of a temperature inversion in the atmospheric boundary layer above the ground. Meteorologically, these freezes are commonly associated with high pressure systems (Winkler et al., 2012b). Terrain influences the near-surface temperatures during radiation freeze events because cold, dense air near the surface flows downhill due to gravity and pools at 9  relatively low elevations (Perry, 1998). In contrast, advection freezes tend to occur with a well-mixed boundary layer and strong winds, and clouds may be present. Advection freezes are commonly associated with cold fronts or the leading edges of polar-origin high pressure systems (Winkler et al., 2012b). Both radiation and advection freezes can significantly decrease fruit production. The type of freeze is critical in preventing cold injury or damage. Radiation type freezes in temperate climates during the spring season tend to be more common than advection freezes (Charrier et al., 2015; Lindkvist et al., 2000; Logan et al., 2000). Since the microclimate influences the temperature during radiation freezes, site selection is a prominent way of preventing freeze damage (Gombos et al. 2011; Perry, 1998), with ideal sites generally on hilltops or slopes (Barden and Neilsen, 2003). Winkler et al. (2012b) examined a radiation freeze event that occurred in 2002 and determined that sour cherry trees located in depressions had relatively higher rates of freeze damage than the trees on hilltops. In contrast, advection freezes may enhance the heat loss and cooling of plant tissue resulting in relatively greater damage, and are more difficult to protect against than radiation type freezes (Perry, 1998). The Winkler et al. (2012b) study concluded that an advection freeze in 2002 caused more damage to sour cherry than those located in relatively poor orchard sites in depressions or other low-lying areas due to greater wind exposure and accompanying sub-freezing temperatures.  A number of studies have compared the frequency and magnitude of freezing temperatures during radiation and advection type events. However, relatively less is known about the relationship between the different freeze types and variations in bud tissue damage. Logan et al. (2000) compared radiation and advection freezes at a peach orchard in Tennessee over two 10  growing seasons. They determined that while radiation freezes occurred more frequently, advection freezes were colder and lasted for longer periods. Similarly, Lindkvist et al. (2000) found that radiation freezes accounted for over 90% of the freezes during the growing season in a Scandinavian mountain range. They concluded that local topography was very important due to the effects of cold air drainage during radiation type freezes. Additionally, a study in the Carpathian Basin in Hungary revealed that during calm, clear nights (i.e. radiation type events) the temperature in a valley setting averaged 3-5°C colder than on nearby slopes (Gombos et al., 2011). Furthermore, Madelin and Beltrando (2005) used topographic features to map frost risk under calm, clear conditions in a vineyard in the Champagne region of France, and the resulting map represented frost damage well except in areas where phenology was delayed due to differing cultivars or varieties. These studies underscore the importance of microclimate and site selection in reducing the risk of freeze damage in commercial fruit production systems. 2.1.3 Influence of Phenology on Freeze Damage Another important factor in determining the risk of springtime freeze damage is the rate of development and phenological stage of the crop. Because fruit trees become increasingly vulnerable to freeze damage as they develop, the stage of bud development is critical in determining the susceptibility a tree (Charrier et al., 2015; Longstroth, 2005). In early developmental stages, temperatures near -6°C begin to cause damage to fruit crops, but as the tree approaches bloom, temperatures as warm as -2°C can cause considerable damage (Longstroth, 2005). However, the role of phenology can be difficult to quantify as the rate of growth and development may vary by cultivar (Chmielewski, 2013; Moghadam et al., 2009; San Martino et al., 2005). Also, because phenology is strongly dependent on temperature, 11  microclimate can affect growth as cooler temperatures may delay development (Jackson, 1966; Logan et al., 2000).  Crop models are frequently used to investigate the complexities of phenology. There are two main associated numerical approaches to simulating phenology: empirical (i.e. statistical) and deterministic (i.e. mechanistic). Empirical-based phenology models correlate timing of plant phenology to climatic factors (Chuine, et al., 2013). The model parameters are estimated from empirical data using various statistical fitting methods, commonly using heat units based on known plant base temperature thresholds (Zhao et al., 2013). Empirical-based phenology models may have limited application and predictive abilities, especially under the context of climate change (Zhao et al., 2013). Deterministic phenology models describe the cause-effect relationships between biological processes and driving environmental factors (Chuine, et al., 2013; Zhao et al., 2013). Parameters of deterministic models in principle can be directly measured; however, this is rarely possible and thus many parameters are estimated using experimental approaches and statistical model-fitting techniques (Chuine, et al., 2013). The deterministic models have more realistic assumptions and can be used in broader applications than empirical models; however, the understanding of the phenological mechanisms may be limited (Zhao et al., 2013). 2.1.4 Frost Risk of Perennial Crops Because springtime freezes can significantly limit temperate perennial fruit production, a number of studies evaluate springtime freeze risk for perennial trees around the world. Most of these studies estimate freeze damage risk by comparing the date of a certain phenological stage to the last date of frost in the spring season (e.g. Fitchett et al., 2014; Kurlus et al., 2013; Ladányi et al., 2010; Molitor et al., 2014). Kurlus et al. (2013) examined the spring frost risk of sour 12  cherry trees in Poland by comparing the date of the last frost to the date of the beginning of flowering. Frost risk was present in most years as the last frost of the spring occurred before the beginning of bloom in only 4 out of the 26 years examined. Furthermore, for the majority of the years, frost occurred within 10 days of the beginning of bloom.  A similar but more complex study by Fitchett et al. (2014) examined frost risk for citrus fruits in two locations in Iran. In addition to last seasonal frost dates and dates of peak flowering, they examined the number and magnitude of frost events occurring after the date of peak flowering. The results at the two locations varied, with large increases in frequency of frost events over time at one location, and decreases over time at the other location. They concluded that the date of peak flowering was advancing more quickly than the last frost date, resulting in an overall increase of frost damage risk. This study also examined the magnitude of temperature during frosts and determined that minimum temperatures had increased over time, reducing the severity of frost damage. Likewise, Ladányi et al. (2010) examined the frequency of frost days and the magnitude of minimum temperatures during a 10 day period before and after bloom to examine the potential damage to sour cherries in Hungary. An average of 1.4 frost days occurred per year, and the number of frost days decreased during the period 1984-2005.  All of the previously mentioned studies considered the impacts of freezes utilizing observations of phenological growth and development. As described above, phenological stages can also be estimated with representative environmental data. For example, Molitor et al. (2014) estimated the date of budburst in the Luxembourgish winegrowing region and compared it to the date of the last frost. They simulated the break of dormancy using a thermal forcing model which accounted for photoperiod and chilling hours, and used a set period of 60 days after the date of 13  budburst to estimate the risk of frost damage across multiple stations. Eleven out of 30 seasons had at least 1 frost event during 1961-1990, and up to 8 frost events occurred during a season. Other studies used more complex methodologies to examine freeze damage by simulating damage for each phenological stage via critical stage-dependent temperature thresholds (Cittadini et al., 2006; Öztekin, 2008). However, the number of studies employing this method is limited due to the relatively small number of field trials or growth chamber studies conducted for specific varieties or cultivars (Winkler et al., 2013). For example, Cittadini et al. (2006) developed a method for examining frost damage risk in sweet cherry trees in South Patagonia in South America. Phenology for multiple cultivars of sweet cherry trees was simulated using cumulative growing degree days (GDDs) beginning on July 15. To estimate frost damage, lethal temperature thresholds for five various phenological stages were used, but these thresholds were not validated for the cultivars or locations in South Patagonia. The risk of frost damage was calculated as the number of seasons during which at least one lethal frost occurred, divided by the 50 years studied, for each location. The results varied from 0.16 to 0.48 depending on location, but there was little difference in frost risk between cultivars. However, the resulting frost damage was not validated.  Similarly, Öztekin (2008) analyzed frost damage in peach trees in Tokat, Turkey, and the results were validated with observations. A phenology model was used in combination with observations of the bloom date or fruit set date each year. Temperature thresholds from a chamber study were used to estimate crop damage amounts of 10%, 50%, and 90%. Crop damage less than 10% was assumed negligible. Greater bud damage was estimated by linear interpolation between the 50% and 90% damage threshold temperatures. They calculated the frequency of damage events by phenological stage, and used the resulting damage amounts to 14  estimate cumulative yield losses for each year. Frost damage occurred in approximately one third of the years, and most of the years with damage had total damage amounts greater or equal to 57%. Almost all of the damage events occurred in the early morning hours, with radiation type freezes. Good agreement was found between estimated and observed yields. Other than the studies mentioned above, quantification of springtime freeze damage on fruit crops is rare in the literature. Additionally, little is known about the physical characteristics of the freeze events themselves. Furthermore, the literature is geographically diffuse, with little information about damaging springtime freezes in the U.S.  2.1.5 Study Objectives The overall goal of this study is to investigate the historical relationship between springtime freezes and sour cherry production in Michigan during the period 1960-2015. The study examines the frequency and spatial variability of springtime freeze events, and assesses the timing of damaging springtime freeze events with respect to sour cherry phenology using a sour cherry simulation model. Finally, the influence of the type and spatial density of various climate-related input data on simulated sour cherry yields is explored using several point-based and gridded climatic datasets.    15  2.2 Data and Methods Weather data used in this analysis of the impact of springtime freeze events on sour cherries were taken from a variety of sources for the period 1960-2015. Observing sites from these sources are shown in Figure 2.1. Observing sites were chosen to be geographically representative of major cherry production areas but also included inland sites further away from Lake Michigan to examine the influence of lake proximity. Hourly temperature, relative humidity, cloud cover, and wind speed data and daily snow cover data were obtained from the NOAA Local Climatological Data (LDD) series for three sites: Traverse City in northwest Lower MI, Muskegon in west central Lower MI, and South Bend in northwest IN. Daily maximum and minimum temperature for the 3 previously mentioned sites and 20 additional surrounding sites from NWS Cooperative Observer Program (COOP) (NOAA Summary of the Day data series) were used to examine the spatial distribution of temperatures during freezes, with 7 sites in the northwest, 8 in the west central, and 8 in the southwest regions. Daily data from 17 NWS COOP sites for 1960-2015 were used in the sour cherry model, with 3 sites, Maple City, Hart, and Eau Claire, used to represent cherry production regions. Additional weather observations were obtained for 2005-2015 from 12 automated stations, which are part of the Michigan Enviroweather Mesonet (www.enviroweather.msu.edu) through Michigan State University Extension, to use in the sour cherry simulation model.  2.2.1 Springtime Freeze Types Spring (March, April, and May) freeze events during the study period were classified at each hourly site as either radiative or non-radiative in nature. Radiative freeze events were identified as freeze events which occurred during the hours 6:00-12:00 UTC (1-7 AM local time) when the observed air temperature reached 0ºC or less, average wind speed was less than or   16   Figure 2.1: Study Weather Station Sites  Weather data observation sites used in the study. 20 NWS COOP stations (1960-2015) with daily data used to examine spatial and temporal variability in springtime freeze events are shown in purple circles, with the purple triangles and associated labels showing the stations with hourly data. Additional stations used in the sour cherry model, with 17 from NWS COOP (1960-2015) and 12 from the Michigan Enviroweather Mesonet (2005-2015), are shown in pink circles. The NWS COOP stations that were used in both analyses are shown in dark purple circles. The 3 reference COOP stations (1960-2015) used in the sour cherry model are labelled and shown in black stars. The red areas represent cherry production acreage from NASS CropScape.  17  equal to 1.54 m/s (3 knots), and sky cover was clear or with few or scattered clouds, with the allowance of one hour with greater cloud cover. This classification scheme is similar to the definition used in Perry (1998). Due to the limited number of hourly data, it was assumed that when a radiative freeze or non-radiative freeze occurred at an hourly station, the same freeze type occurred at the surrounding nearby COOP sites. The average minimum temperatures during freeze events were interpolated using an inverse distance weighting scheme and mapped to identify any spatial trends over time. Days with missing data during the 6:00-12:00 UTC time interval were omitted. 2.2.2 Sour Cherry Simulation Model  A sour cherry simulation model developed by Black et al. (in preparation) was utilized to examine the relationship between sour cherry crop damage and springtime freeze events. The model is a hybrid empirical/process-based model that requires input of daily maximum and minimum temperatures. The model is based on the work of Zavalloni et al. (2006) to simulate sour cherry phenology using seasonally summed base 4ºC GDDs. Seasonal GDD accumulations began on March 1 of each year, and eight different phenological stages were simulated (Table 2.1). Critical temperature thresholds for each phenological stage based on a chamber study by Dennis and Howell (1974) were used to estimate damage due to freezes. In the model, the freeze damage during each phenological stage was accumulated throughout the growing season, resulting in a summed viable buds remaining index. Buds remaining represents the proportion of the buds that survived springtime freeze damage, with a value of 0 indicating that no buds survived, and a value of 1 indicating all buds survived. The simulated phenology used in the sour cherry model was validated in Michigan by Zavalloni et al. (2006). The sour cherry model output of buds remaining at Maple City was validated with annual observed sour cherry yield data  18   aggregated over five counties for the northwest and west central production regions in Lower Michigan for 1982-2015 from United States Department of Agriculture (USDA) National Agricultural Statistics Service (NASS). The results were satisfactory with the model accounting for 30% of the variation in the northwest observed yield, and a linear regression between the observed yields and buds remaining was statistically significant at an alpha level of 0.001. In order to better understand the impact of springtime freezes on sour cherry production on bud cold damage, a simulation model was used with historical climatic data to examine temporal patterns and trends. Maple City, Hart, and Eau Claire, MI, representing Northwest, West Central, and Southwestern production areas, respectively were selected as representative regional sites in this portion of the study. Daily climatic data required by the sour cherry model were obtained from the NOAA Summary of the Day series for the period 1960-2015.  2.2.3 Sour Cherry Simulations over Space  Sour cherry yields are known to vary widely over both time and space, largely the result of microclimatic differences over relatively short distances (Winkler et al., 2013). In an attempt to examine spatial yield variability over time, freeze impacts were examined in the  Table 2.1: Base 4ºC GDD Accumulations for Phenological Stages GDD accumulation totals used in the sour cherry simulation model to define the phenological stages of the crop, with the GDD thresholds adapted from Zavalloni et al. (2006). Phenological Stage GDD4 Accumulation from March 1 0 0-120 2 Side Green 120-153 3 Green Tip 154-173 4 Tight Cluster 174-189 5 Open Cluster 190-207 6 First White 208-227 7 First Bloom 228-252 8 Full Bloom 253-310 9 Petal Fall 311+  19  northwest and west central Lower Michigan production areas using various types of input climatic data, described below, in the sour cherry simulation model (Table 2.2). The output of buds remaining was interpolated over the regions using kriging. Cherry acreage data obtained from the NASS CropScape 30 m resolution cropland data layer (USDA NASS, 2013) was used to spatially weight the interpolated output. (See Appendix for additional information on NASS CropScape.) Production area data from the 2013 growing season was used for all model runs, as the areas of cherry production have not changed significantly in the past few decades. The spatially weighted buds remaining output was then summed for the northwest and west central regions, and the resulting time series were compared to the observed regional yields for northwest and west central Lower Michigan to determine if the use of such a weighted input data set in the sour cherry model could provide more accurate estimations of yield over a region, versus from a single representative point. The NWS COOP stations Maple City and Hart were used in the sour cherry model to examine simulated model output of buds remaining for 1960-2015 at individual locations. These int data were obtained from weather observations for 2005-2015 from 29 stations across Michigan, using data from 12 automated stations which are part of the Michigan Enviroweather Mesonet (www.enviroweather.msu.edu) through Michigan State University Extension and 17 station sites from the NWS COOP network as it theoretically provides additional information into the regional estimation process.  To account for spatial variation over a longer time period, a hybrid gridded dataset was developed using daily observations from a single site, Maple City and from the additional site data from the Michigan Enviroweather Mesonet and NWS COOP for 1960-2015. The data series 20  Table 2.2: Various Climatic Input Data for the Sour Cherry Simulation Model Various climatic input data types used in the sour cherry simulation model and their date ranges. Name Input Data Type Description Range Maple City Data Single Point Single point data from the NWS COOP station Maple City 1960-2015 Hart Data Single Point Single point data from the NWS COOP station Hart 1960-2015 Additional Station Data Multiple Point Multiple point data from 17 NWS COOP and 12 Enviroweather Mesonet sites 2005-2015 Side Green Data Hybrid Gridded Single point data from Maple City converted into a gridded surface using multiple point data, corrected by phenological stage side green 1960-2015 Bloom Data Hybrid Gridded Single point data from Maple City converted into a gridded surface using multiple point data, corrected by phenological stage bloom 1960-2015 PRISM Data Gridded Gridded data with a 4 km spatial resolution 1981-2014 NLDAS Data Gridded Gridded data with a 1/8° (12.5 km) spatial resolution 1979-2015   at Maple City was temporally adjusted to represent the weather at other locations in an attempt to reflect known sub-regional differences in phenology. The average differences in simulated phenology dates were calculated between Maple City and the other locations for the key phenological stages, side green (stage 2) and full bloom (stage 8). The dates of observations at Maple City were then shifted corresponding to the average differences in phenological dates for 2005-2015 between Maple City and the other locations. This process resulted in two new datasets for each of the 28 station sites for 1960-2015, with one dataset based on differences in he other based on differences in dates .   Existing gridded climatic data were also used to attempt to better capture spatial variability in sour cherry yields. Daily maximum and minimum temperature data for the 29 21  weather station locations were extracted from two gridded datasets, the Parameter-elevation Regressions on Independent Slopes Model (PRISM) data set provided by the PRISM Climate Group at Oregon State University (Daly et al., 2008; PRISM Climate Group) and the North American Land Data Assimilation System Phase 2 data set (NLDAS-2) (Mitchell et al., 2004). The PRISM dataset was developed through the use of a weighted linear regression to interpolate point climate data as a function of elevation across the U.S. (Daly et al., 2008). The PRISM dataset accounts for the major physiographic factors influencing climate patterns through a digital elevation model, including location, elevation, topographic orientation, and cold air drainage (Daly et al., 2008). The data are available at a daily temporal resolution for 1981-2014 and have a 4 km spatial resolution.   Unlike the PRISM dataset, the NLDAS-2 dataset does not rely on physiographic factors, and instead integrates observation based data with model reanalysis data (Mitchell et al., 2004). NLDAS-2 relies on the North American Regional Reanalysis (NARR) dataset, mainly for non-precipitation fields, and temperatures are adjusted with satellite-derived data (Xia et al., 2012). The precipitation amounts are dominated by observed data, and when possible, the observations are interpolated using the PRISM technique and temporally disaggregated using radar data (Xia et al., 2012). NLDAS-2 data have an hourly temporal resolution, 1/8° (~12.5 km) spatial resolution, and are available for 1979-2015.  In this study, the nearest grid cell values for station locations were extracted from the gridded climate datasets, and the grid cells containing a majority of water were omitted from the analysis. The use of the .   22  2.3 Results and Discussion The frequency and spatial variation of springtime freezes, their physical characteristics, and timing with respect to crop phenology are important for understanding and predicting sour cherry production. Additionally, information on radiative versus non-radiative freezes can help growers manage their crops efficiency as freeze protection methods depends on the type of freeze (Perry, 1998). 2.3.1 Temporal and Spatial Variability of Springtime Freezes The frequencies of freeze types were classified using hourly data at Traverse City, Muskegon, and South Bend; however, due to the lack of hourly data in the region, it was assumed that when a freeze occurred at one of the hourly stations, the same freeze type occurred at the COOP stations with daily data in the associated region. However, this is a potential source of error in the analysis as the regions are not completely homogeneous over space. To better quantify this potential error, the number of days when a freeze was classified at an hourly station and the minimum temperature at a COOP station did not reach freezing were summed. The results varied across the region, indicating some lack of spatial consistency. The majority of stations had less than 5% of days with above freezing minimum temperatures during radiation or non-radiation classified springtime freeze events. However, a few stations had larger totals, with up to 29% of regional radiation freezes (168 out of 573 days) not observed at Frankfort, which is located in the northwestern Lower Peninsula on the shore of Lake Michigan. These results highlight the limitation of choosing individual stations to represent regions, and the uncertainty in the calculating average minimum temperatures during freeze events. The frequency of all springtime freeze events and radiation events were examined for the three hourly stations Traverse City, Muskegon, and South Bend for 1960-2015 (Table 2.3). Only 23  the months of March, April, and May were analyzed, as no freeze events occurred in June or later summer months. As expected, the frequency of springtime freeze events was greatest at Traverse City, the northernmost location, with a total of 2,015 events over the period of record, followed by Muskegon with 1,596, and then South Bend with 1,111 events. The majority of springtime freeze events occurred in March, followed by April and then May at all locations. These results are consistent with existing climatic normals, latitudinal gradients, and the annual temperature cycle.  Only a small proportion of the springtime freezes, ranging from 9 to 35% by location, were classified as radiation type events in March and April, according to our classification criteria of calm, clear atmospheric conditions. In May, the proportion of radiation freezes was slightly larger, with 55% at Traverse City and approximately 40% at Muskegon and South Bend. This result differs somewhat from previous studies conducted in Tennessee, U.S. and Sweden which indicated that radiation freezes were the most common freeze type (Lindkvist et al., 2000; Logan et al., 2000). The difference in results may be accounted for by the differing classification scheme and much longer period of record used in this study. The data also shows that the frequency of radiation freezes (and the proportion of all freezes that are the radiation type) Table 2.3: Frequency of Freezes at Traverse City, Muskegon, and South Bend, 1960-2015 Frequency of freeze events for spring season and by month at Traverse City, Muskegon, and South Bend, 1960-2015.   All Freezes (n) Radiation Freezes (n)   Traverse City Muskegon South Bend Traverse City Muskegon South Bend March 1152 1088 809 237 100 111 April 655 448 287 227 88 79 May 208 60 15 114 25 6 Total 2015 1596 1111 578 213 196  24  increased monthly from March to May at all locations. Finally, springtime radiation freeze events occurred simultaneously at the three sites on only 41 days between 1960 and 2015, highlighting the considerable spatial variability over time across the region. The radiation freeze events for each region are classified further as those with or without more than 2.5 cm of snow cover (Table 2.4). Snow cover is known as an effective insulator of the ground surface and is typically associated with relatively colder minimum surface temperatures. At all locations the majority of radiation freeze events with snow cover occurred in March with none occurring in May. Traverse City had the highest percentages of radiation freezes with snow cover during spring at 30%, followed by Muskegon at 19%, and then South Bend at 8%. The percentage of spring radiation freezes with snow cover decreased over time during the season at all stations. At Traverse City and Muskegon, 64% and 35% of the radiation freezes were associated with snow cover, respectively, but in April only 10% and 6% were, respectively. In South Bend, 14% of the radiation freeze events in March had snow cover, with no radiation freeze events with snow cover in April. Table 2.4: Frequency of Freezes with 2.5+ cm of Snow Cover at Traverse City, Muskegon, and South Bend, 1960-2015  Frequency of freezes with at least 2.5cm of snow cover for each month and spring season total at Traverse City, Muskegon, and South Bend, 1960-2015.  Freezes with Snow Cover (n) Radiation Freezes with Snow Cover (n)  Traverse City Muskegon South Bend Traverse City Muskegon South Bend March 237 100 111 151 35 15 April 227 88 79 22 5 0 May 114 25 6 0 0 0 Total 578 213 196 173 40 15   25  The potential impacts of the freeze events on cherry cold injury at each reference site were examined with the sour cherry simulation model (Table 2.5). A damaging freeze event was defined as a day when the model simulated a loss of buds due to cold temperature. As expected, Traverse City had the most damaging freezes, with 100 events during spring over the 56 years, followed by South Bend with 38 events, and then Muskegon with 23. By month, the majority of damaging freezes at Traverse City and Muskegon occurred in May, but at South Bend the most frequent month was April. In March no damaging freezes occurred at Traverse City or Muskegon, and only 1 occurred at South Bend. This pattern is a result of crop phenology, as the risk of freeze damage increases as the crop advances through bloom. The majority of the damaging freezes were radiation freeze types at Traverse City (67%) and at Muskegon (61%), but only 34% of the damaging events at South Bend were radiation freezes. Because microclimate plays a significant role on temperature during freeze events, the average minimum temperatures during the springtime freeze events for 1960-2015 were examined across western Michigan and northwest Indiana. The average minimum temperatures for radiation and non-radiation freezes are shown by month for each station in Table 2.6. The  Table 2.5: Frequency of Simulated Damaging Freezes at Traverse City, Muskegon, and South Bend, 1960-2015 Frequency of simulated damaging freezes and freeze types at Traverse City, Muskegon, and South Bend by month and for the spring season, 1960-2015.   Simulated Damaging Freezes (n) Simulated Damaging Radiation Freezes (n)   Traverse City Muskegon South Bend Traverse City Muskegon South Bend March 0 0 1 0 0 1 April 36 8 35 25 5 11 May 64 15 2 42 9 1 Total 100 23 38 67 14 13  26  Table 2.6: Average Minimum Temperature during Radiation and Non-Radiation Freeze Events, 1960-2015 Average minimum temperatures (°C) using daily data from the NWS COOP stations for the spring months, during radiation and non-radiation freeze events classified from hourly data for 1960-2015.     March April May Seasonal     Rad Non-Rad Rad Non-Rad Rad Non-Rad Rad Non-Rad Region Station Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Northwest Cadillac -12.5 231 -9.8 884 -4.8 221 -4.6 420 -2.2 112 -2.3 94 -7.4 564 -7.8 1398 East Jordan -11.8 236 -9.0 909 -4.2 225 -4.3 425 -2.1 114 -2.0 94 -6.9 575 -7.1 1428 Frankfort -6.5 235 -6.1 904 -0.6 226 -2.2 418 1.3 112 0.3 92 -2.6 573 -4.6 1414 Ludington -7.6 203 -6.1 834 -2.3 209 -2.4 402 0.0 106 -0.5 93 -3.9 518 -4.6 1329 Manistee -8.0 234 -6.0 888 -1.7 218 -1.8 392 0.5 111 0.6 92 -3.9 563 -4.3 1372 Maple City -10.3 237 -7.6 915 -4.2 227 -3.7 428 -2.5 114 -2.2 94 -6.4 578 -6.1 1437 Traverse City -11.0 237 -7.9 915 -4.3 227 -3.6 428 -2.3 114 -1.9 94 -6.7 578 -6.2 1437 Regional Average -9.7 230.4 -7.5 892.7 -3.2 221.9 -3.3 416.1 -1.1 111.9 -1.1 93.3 -5.4 564.1 -5.8 1402.1 West Central Big Rapids -10.1 97 -9.2 966 -4.9 87 -4.9 352 -2.0 25 -2.5 34 -7.0 209 -7.9 1352 Grand Rapids -7.8 100 -6.4 988 -3.2 88 -3.0 360 -0.6 25 -1.5 35 -5.0 213 -5.4 1383 Greenville -9.7 100 -7.9 987 -4.3 88 -4.4 360 -2.0 25 -2.7 35 -6.6 213 -6.9 1382 Hart -9.0 100 -7.2 988 -4.0 88 -3.3 360 -2.2 25 -1.9 35 -6.1 213 -6.0 1383 Hastings -9.1 99 -7.2 976 -4.2 84 -3.7 345 -1.7 24 -1.9 35 -6.3 207 -6.1 1356 Holland -7.9 89 -5.7 925 -3.4 76 -2.7 346 -1.2 25 -1.7 33 -5.2 190 -4.8 1304 Montague -10.1 97 -6.9 926 -6.3 85 -4.4 339 -3.5 23 -3.2 33 -7.8 205 -6.1 1298 Muskegon -8.5 100 -6.0 988 -4.0 88 -2.7 360 -2.0 25 -1.3 35 -5.9 213 -5.0 1383 Regional Average -9.0 97.8 -7.1 968.0 -4.3 85.5 -3.6 352.8 -1.9 24.6 -2.1 34.4 -6.2 207.9 -6.0 1355.1    27   Southwest Benton Harbor -8.1 109 -5.8 678 -4.3 75 -3.2 196 -2.3 6 -1.5 9 -6.4 190 -5.2 883 Dowagiac -8.1 107 -6.9 685 -4.2 78 -4.1 207 -0.7 6 -3.0 9 -6.3 191 -6.2 901 Goshen -6.1 111 -5.4 689 -2.2 76 -2.7 197 -0.3 6 -1.3 9 -4.4 193 -4.7 895 Lagrange -7.3 96 -6.7 607 -3.4 67 -3.5 168 -0.2 6 -2.5 9 -5.5 169 -5.9 784 La Porte -4.7 107 -4.9 653 -0.5 76 -1.7 197 1.9 6 0.3 9 -2.8 189 -4.1 859 South Bend -6.7 111 -5.4 698 -2.8 79 -2.9 208 -1.3 6 -1.4 9 -4.9 196 -4.8 915 Three Rivers -7.5 111 -6.7 688 -3.2 77 -3.6 207 -1.4 6 -2.5 9 -5.6 194 -5.9 904 Wantanah -7.2 106 -6.2 672 -4.1 72 -3.1 195 -1.5 6 -1.9 9 -5.8 184 -5.4 876 Regional Average -7.0 107.3 -6.0 671.3 -3.1 75.0 -3.1 196.9 -0.7 6.0 -1.7 9.0 -5.2 188.3 -5.3 877.1  springtime average minimum temperatures during radiation freezes range from -7.8 to -2.6°C across stations while the non-radiation freezes were slightly cooler, ranging from -7.9 to -4.1°C. Lack of spatial consistency across regions is shown, and some stations are consistently colder or warmer than others. For example, Frankfort, in the close vicinity of Lake Michigan, is warmer on average than the surrounding stations during radiation freezes. The radiation freezes are generally colder than the non-radiation freezes in March, with temperatures more similar between freeze types in April and May. The average minimum temperatures for all freeze events, radiation freezes, and non-radiation freezes in March are shown in Figure 2.2. The mean temperatures are relatively warmer in the southwest and colder in the northeast, with a    28  change from -9.3°C to -4.4°C. This temperature gradient is likely due to the combined influences of latitude (a regional difference of 3.7°), the lake influence, topography, and the various microclimates of the individual sites. This result highlights the importance of choosing a suitable orchard site in order to reduce the overall risk of springtime freezes. On average, the temperatures are 1-2°C colder during radiation freezes than non-radiation freezes in March across all regions. This result differs from the Logan et al. (2000) study which determined that advection freezes were colder than radiation freezes. Again this may be due to the different location or classification criteria, as this study categorized non-radiation freezes instead of advection freezes (some of which could have been radiation-like in nature).  The average minimum temperatures during radiation freeze events with and without snow cover are shown by station in Table 2.7. The average minimum temperatures of springtime radiation freezes range from -14.3 to -7.8°C for those with snow cover and -6.4 to -0.5 °C without snow cover. In general, the springtime minimum temperatures during radiation freezes are 7-8°C colder with snow cover than without across the three regions. Again, the temperatures vary by station, even within region, up to a difference in 6.5°C during radiation freezes with snow cover in the northwest, highlighting microclimatic effects. The effects of microclimate were further highlighted in the interpolated map of the temperatures (Figure 2.3). Microclimate influences were especially apparent in the northwest region, as the average minimum temperatures during radiation freezes in March were warmer near Lake Michigan than the locations further inland. The radiation freeze events with snow cover were 5-6°C colder on average in March than the events without snow cover across all regions.  The average minimum temperatures during simulated damaging radiation and non-radiation freeze events are shown in Table 2.8. These averages are based on smaller numbers of 29  freezes, and again the temperatures vary by station. The average minimum temperatures range from -5.5 to 0.3°C during radiation freeze types and from -6.4 to -0.4°C during non-radiation freeze types. The differences in average minimum temperatures between damaging radiation and non-radiation freeze events are small, and are not consistent across month or region. 30   Figure 2.2: Interpolated Average Minimum Temperature of March Freeze Events, 1960-2015  Average minimum temperature (°C) for of all March freeze events (left panel) and for radiation (center) and non-radiation (right) type freezes, 1960-2015.  All freezes included all days with freezing temperatures at each station with daily data, and the days included in the radiation and non-radiation freezes were classified based on hourly data at the three reference stations.  31    Table 2.7: Average Minimum Temperatures during Radiation Freeze Events with 2.5+ cm of Snow Cover, 1960-2015 Average minimum temperatures (°C) using daily data from the NWS COOP stations for the spring months, during radiation freeze events with and without 2.5+ cm of snow cover classified from hourly data for 1960-2015.     March April May* Seasonal     With Snow Without Snow With Snow Without Snow Without Snow With Snow Without Snow Region Station Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Northwest Cadillac -14.9 151 -8.1 80 -10.2 22 -4.2 199 -2.2 112 -14.3 173 -4.4 391 East Jordan -14.8 150 -6.5 86 -9.7 22 -3.6 203 -2.1 114 -14.1 172 -3.8 403 Frankfort -8.2 149 -3.4 86 -4.7 22 -0.2 204 1.3 112 -7.8 171 -0.5 402 Ludington -9.9 120 -4.4 83 -7.3 16 -1.9 193 0.0 106 -9.6 136 -1.9 382 Manistee -9.8 149 -4.9 85 -5.4 19 -1.3 199 0.5 111 -9.3 168 -1.6 395 Maple City -12.4 151 -6.7 86 -9.0 22 -3.7 205 -2.5 114 -11.9 173 -4.0 405 Traverse City -13.4 151 -6.9 86 -9.2 22 -3.8 205 -2.3 114 -12.8 173 -4.0 405 Regional Average -12.0 145.9 -5.8 84.6 -8.0 20.7 -2.7 201.1 -1.1 111.9 -11.5 166.6 -2.9 397.6 West Central Big Rapids -13.6 34 -8.2 63 -10.9 5 -4.6 82 -2.0 25 -13.3 39 -5.5 170 Grand Rapids -10.9 35 -6.1 65 -10.2 5 -2.7 83 -0.6 25 -10.8 40 -3.7 173 Greenville -13.6 35 -7.7 65 -10.2 5 -4.0 83 -2.0 25 -13.2 40 -5.1 173 Hart -13.3 35 -6.7 65 -11.8 5 -3.5 83 -2.2 25 -13.1 40 -4.5 173 Hastings -12.3 35 -7.3 64 -11.0 5 -3.8 79 -1.7 24 -12.1 40 -4.8 167 Holland -11.1 28 -6.4 61 -10.1 5 -3.0 71 -1.2 25 -10.9 33 -4.0 157 Montague -13.3 34 -8.4 63 -14.8 5 -5.7 80 -3.5 23 -13.5 39 -6.4 166 Muskegon -11.8 35 -6.7 65 -12.3 5 -3.5 83 -2.0 25 -11.9 40 -4.5 173 Regional Average -12.5 33.9 -7.2 63.9 -11.4 5.0 -3.9 80.5 -1.9 24.6 -12.4 38.9 -4.8 169.0      32  Table 2.7  Southwest Benton Harbor -12.3 14 -7.5 95 - 0 -4.3 75 -2.3 6 -12.3 14 -6.0 176 Dowagiac -13.2 14 -7.3 93 - 0 -4.2 78 -0.7 6 -13.2 14 -5.7 177 Goshen -11.0 15 -5.3 96 - 0 -2.2 76 -0.3 6 -11.0 15 -3.8 178 Lagrange -13.2 12 -6.5 84 - 0 -3.4 67 -0.2 6 -13.2 12 -4.9 157 La Porte -9.3 14 -4.0 93 - 0 -0.5 76 1.9 6 -9.3 14 -2.3 175 South Bend -11.8 15 -5.9 96 - 0 -2.8 79 -1.3 6 -11.8 15 -4.4 181 Three Rivers -12.5 15 -6.7 96 - 0 -3.2 77 -1.4 6 -12.5 15 -5.0 179 Wantanah -12.9 15 -6.3 91 - 0 -4.1 72 -1.5 6 -12.9 15 -5.2 169 Regional Average -12.0 14.3 -6.2 93.0 - 0.0 -3.1 75.0 -0.7 6.0 -12.0 14.3 -4.6 174.0 *No radiation freeze events with snow cover in May   33   Figure 2.3: Interpolated Average Minimum Temperature of March Radiation Type Freeze Events, 1960-2015  Average minimum temperature (ºC) of March radiation type freeze events (left panel), radiation type freezes with snow cover (center), radiation type freezes without snow cover (right), 1960-2015. The days included were classified based on hourly data at the three reference stations.    34  Table 2.8: Average Minimum Temperatures during Simulated Damaging Freeze Events, 1960-2015 Average minimum temperatures (°C) using daily data from the NWS COOP stations for the spring months, during radiation and non-radiation simulated damaging freeze events through the sour cherry model for 1960-2015.     March* April May Seasonal     Radiation Radiation Non-Radiation Radiation Non-Radiation Radiation Non-Radiation Region Station Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Tmin (°C) (N) Northwest Cadillac - 0 -4.4 25 -4.6 11 -3.3 40 -2.8 25 -3.7 65 -3.3 36 East Jordan - 0 -3.7 25 -3.9 10 -3.3 42 -3.0 25 -3.5 67 -3.2 35 Frankfort - 0 0.1 25 -2.0 11 0.4 42 0.1 25 0.3 67 -0.6 36 Ludington - 0 -2.9 25 -1.8 11 -1.0 40 -1.3 25 -1.7 65 -1.5 36 Manistee - 0 -1.9 25 -1.5 11 -0.1 41 0.2 25 -0.8 66 -0.4 36 Maple City - 0 -4.1 25 -4.0 11 -3.6 42 -2.8 25 -3.8 67 -3.2 36 Traverse City - 0 -3.5 25 -4.3 11 -3.4 42 -3.2 25 -3.5 67 -3.6 36 Regional Average - 0.0 -2.9 25.0 -3.2 10.9 -2.0 41.3 -1.8 25.0 -2.4 66.3 -2.2 35.9 West Central Big Rapids - 0 -5.6 5 -4.8 3 -3.1 9 -4.7 5 -4.0 14 -4.7 8 Grand Rapids - 0 -3.3 5 -3.2 3 -0.6 9 -2.8 6 -1.6 14 -2.9 9 Greenville - 0 -4.8 5 -5.4 3 -1.7 9 -4.4 6 -2.8 14 -4.7 9 Hart - 0 -3.4 5 -4.3 3 -2.7 9 -4.4 6 -3.0 14 -4.3 9 Hastings - 0 -3.4 5 -4.4 3 -1.5 9 -2.8 6 -2.2 14 -3.3 9 Holland - 0 -1.7 4 -4.4 2 -1.7 9 -3.3 6 -1.7 13 -3.6 8 Montague - 0 -5.2 5 -7.8 3 -3.4 9 -5.6 5 -4.1 14 -6.4 8 Muskegon - 0 -3.8 5 -3.7 3 -2.8 9 -3.0 6 -3.1 14 -3.2 9 Regional Average - 0.0 -4.0 4.9 -4.8 2.9 -2.2 9.0 -3.8 5.8 -2.8 13.9 -4.1 8.6     35  Table 2.8  Southwest Benton Harbor -6.7 1 -5.5 11 -4.2 24 -3.9 1 -0.6 1 -5.5 13 -4.1 25 Dowagiac -5.6 1 -5.2 11 -5.3 24 -2.2 1 -3.9 1 -5.0 13 -5.2 25 Goshen -2.2 1 -2.7 10 -3.8 24 -1.1 1 -2.2 1 -2.6 12 -3.8 25 Lagrange -3.3 1 -3.3 7 -4.6 16 -3.3 1 -2.2 1 -3.3 9 -4.4 17 La Porte -1.7 1 -0.6 11 -2.3 21 -1.1 1 0.0 1 -0.7 13 -2.2 22 South Bend -5.6 1 -3.4 11 -4.1 24 -3.3 1 -2.2 1 -3.6 13 -4.0 25 Three Rivers -3.9 1 -3.6 10 -4.3 24 -3.3 1 -3.9 1 -3.6 12 -4.3 25 Wantanah -6.1 1 -4.3 11 -3.7 24 -3.3 1 -2.2 1 -4.4 13 -3.6 25 Regional Average -4.4 1.0 -3.6 10.3 -4.0 22.6 -2.7 1.0 -2.2 1.0 -3.6 12.3 -3.9 23.6 *One radiation freeze and no non-radiation freezes were damaging in March   36  2.3.2 Characteristics of Damaging Freezes To examine the physical characteristics of the simulated damaging freeze events on a diurnal basis, statistics of temperature, wind speed, sky cover, and dewpoint for all the radiation type freezes were calculated for each hour at the three reference stations with hourly data. The hourly characteristics for damaging freezes in April and May at each station were similar, thus only the characteristics of May at Traverse City, which had the highest number of events, are shown. These statistics are presented in Figure 2.4 by hour in the form of box and whisker plots for all simulated damaging radiation freeze events during May at Traverse City from 1960 to 2015 (42 total). During these freeze events, average temperatures tended to decrease rapidly from 22:00 to 06:00 UTC, associated with decreases in wind speed and sky cover. These characteristics were expected as clear, calm nighttime conditions were used in the criteria to classify radiation freezes. There were also relatively small concurrent decreases in dewpoint likely associated with condensation or sublimation and decreases in atmospheric humidity. During the overnight hours, the sharpest hourly decline in temperature on average occurred between 1:00-2:00 UTC with an average decrease of 2°C.  For comparison, the statistical averages of hourly conditions for simulated damaging non-radiation type freezes in May at Traverse City (22 events) are given in Figure 2.5. The non-radiation damaging freezes had higher variability in temperature, wind, and sky cover, and more outliers for dewpoint, where outliers were defined as data points above and below a distance of 1.5 times the interquartile range from the 1st and 3rd quartiles. Additionally, the diurnal cycle for all variables was less distinct for the non-radiation damaging freezes. There was a similar decrease in clouds and wind at night during the non-radiation damaging freezes because relatively clear, calm conditions may still have played a role in the freeze events, although the 37       Figure 2.4: Hourly Characteristics of Simulated Damaging Radiation Type Freeze Events Statistical box and whisker plots of air temperature (ºC), wind speed (m/s), sky cover category, and dewpoint (ºC) for simulated damaging radiation type freeze events in May for 1960-2015 at Traverse City for hours 13:00 to 12:00 UTC. A sky cover value of '1' indicates clear skies, '2' for few clouds, '3' for scattered clouds, '4' for broken sky cover, and '5' is overcast. Box plots include median value (thick black line), 25th and 75th percentiles (bottom and top of boxes), distance of 1.5 times the interquartile range from the 1st and 3rd quartiles (bottom and top whiskers), and series extremes (open circles) for each hour.  38       Figure  Figure 2.5: Hourly Characteristics of Simulated Damaging Non-radiation Type Freeze Events Statistical box and whisker plots of air temperature (ºC), wind speed (m/s), sky cover category, and dewpoint (ºC) for simulated non-radiation type damaging radiation freeze events in May for 1960-2015 at Traverse City for hours 13:00 to 12:00 UTC. A sky cover value of '1' indicates clear skies, '2' for few clouds, '3' for scattered clouds, '4' for broken sky cover, and '5' is overcast. Box plots include median value (thick black line), 25th and 75th percentiles (bottom and top of boxes), distance of 1.5 times the interquartile range from the 1st and 3rd quartiles (bottom and top whiskers), and series extremes (open circles) for each hour. 39  actual events did not meet the criteria for a radiation type freeze event. On average, during hours 1:00-12:00 UTC of non-radiation damaging freezes, the greatest drop in temperature occurred from 0:00 to 1:00 UTC with an average decrease of 1.4°C. The timing of model simulated damaging freezes with respect to phenological stage was investigated at Maple City, Hart, and Eau Claire for 1960-2015. The average frequencies of damage events per year are shown in the top row of Figure 2.6. The averages were based on 162 damaging freeze events that occurred at Maple City, 68 events at Hart, and 41 at Eau Claire. On average, stage 8 had the highest average frequency of damage events per year at Maple City (0.63), followed by stage 6 (0.48) and then stage 2 (0.39). Similarly, at Hart stage 8 had the highest average frequency of damage events per year (0.29), followed by stages 5 and 6 (0.20) and then stage 3 (0.16). The pattern at Eau Claire was somewhat different, with the highest average frequency in stage 2 (0.20), followed by stage 5 (0.16), and then stages 3 and 8 (0.11). At Eau Claire, there were no damaging freeze events that occurred during phenological stage 4 from 1960 to 2015 due to the relatively short stage length. To account for the varying lengths of each phenological stage, the frequency totals were normalized by the length of the stage (days) to get average frequency per day in each phenological stage (Figure 2.6; middle row). Across all stations, the later phenological stages 5-8 had the highest average frequency of damaging freeze events per day. This result was expected since the later phenological stages are at a higher risk to freeze damage than the earlier stages given relatively higher damage threshold temperatures. At Maple City, stage 6 had the highest normalized average frequency at 0.178. At Hart, stage 5 had the highest average normalized frequency at 0.080, followed closely by stage 6 at 0.075. Similarly, stage 5 had the highest average normalized frequency at Eau Claire at 0.074. 40   Figure 2.6: Average Frequency and Severity of Simulated Damaging Freezes by Phenological Stage, 1960-2015 Average frequency and severity of simulated damaging freezes with respect to phenological stages 2-9 for 1960-2015 for Maple City, Hart, and Eau Claire. Average frequency of damaging events per year are given in the top row, average frequency of damaging events per day normalized by stage length in the middle row, and average severity of damage per damaging event in the bottom row.     41  While frequency of damaging freezes is important, it is also critical to understand the severity of the freezes. Therefore, the average simulated damage amount per damaging freeze event was examined for each phenological stage (Figure 2.6; bottom row). On average, freeze damage was most severe in the earlier phenological stages, 2-3, for all three stations. At Maple City, the most severe average damage per event occurred in stages 2 and 3 with an average damage severity of 0.27. For Hart the highest average damage amount occurred in stage 2 at 0.25, closely followed by stage 3 at 0.24. At Eau Claire stage 2 had an average damage severity of 0.28 and 0.25 for stage 3. The other phenological stages, 4-9, experienced less severe damage, ranging from 0.11 to 0.18 across the stations. This result of most severe damage in the earlier phenological stages, but more frequent damage in the later stages may be an important consideration for growers to understand the risk of freeze damage as the fruit tree develops and to better protect their crops from springtime freeze damage. 2.3.3 Influence of Various Climatic Input Data on Sour Cherry Simulations Because of the lack of representative, lengthy climatic data series in the study area (and especially directly within the production areas themselves), a number of spatially weighted versions of input data were developed and used in the sour cherry model for the northwest and west central production regions. The representativeness of the various input datasets was assessed by comparing the simulated buds remaining to the observed yields for northwest and west central Michigan. Observed yields were not available for southwest Michigan. The ratio of buds remaining represents the proportion of the buds that survived after accounting for damage from springtime freezes, and therefore should be correlated with observed yields. The model simulated buds remaining is referred to as the output for the various input datasets as a shorthand. The spatially weighted input data variants were the Side Green, Bloom, Additional 42  Station, NLDAS, and PRISM data. The model outputs were compared among the different series and with the output from a single point series at the reference stations for the northwest (Maple City) and west central regions (Hart).  The relationships between model simulated buds remaining using various input datasets and the observed yields for the northwest and west central regions are shown in Table 2.9, using linear regressions. For the northwest region, the linear regression indicates all the outputs from various input data are statistically significant at the 0.05 alpha level, except for the output from NLDAS data. Similarly, all the linear regressions except the NLDAS and the Hart series outputs for the west central region were statistically significant.  Table 2.9: Linear Regression for Observed Yields and Model Output Linear regression with model simulated buds remaining from various input data as the dependent variable and observed yields in the northwest (top) and west central (bottom) regions as the independent variable. Northwest Region Data Input for Model R2 P-value Range Maple 0.3 0.001*** 1982-2015 Side Green 0.31 0.001*** 1982-2015 Bloom 0.33 0.000*** 1982-2015 Additional Station 0.46 0.022* 2005-2015 PRISM 0.21 0.007** 1982-2014 NLDAS 0 0.967 1982-2015 West Central Region Hart 0.07 0.117 1982-2015 Side Green 0.22 0.005** 1982-2015 Bloom 0.22 0.005** 1982-2015 Additional Station 0.49 0.017* 2005-2015 PRISM 0.38 0.000*** 1982-2014 NLDAS 0 0.981 1982-2015 Significance: * for 0.05, ** for 0.01, *** for 0.001    The single point data, Maple City data, accounts for 30% of the variation in the observed yields in the northwest; however, the Hart data only accounts for 7% of the observed yields in the west central region. The outputs from the Side Green and Bloom data improve the regression 43  results slightly for the northwest region, and largely for the west central region. For both the northwest and west central regions, the output from the Additional Station data accounts for the highest percentage of variation in observed yields, 46% and 49%, respectively. The output from the PRISM data lowers the R squared value compared to the single point data output in the northwest, but improves the R squared value in the west central region, accounting for 38% of the observed yields.  When plotted in time series form, the model obtained output from Maple City data represents the observed northwest yields fairly well (Figure 2.7). The observed yields show major decreases in yield in 2002 and 2012, which are years with known springtime freeze events. The Maple City output represents the observed yields well in 2012, but the buds remaining output only decreases to approximately 0.5 in 2002, while the observed yield is extremely low. The Maple City output overestimates damage, which can be seen in 1989 and 2010 when simulated buds remaining are near 0 but the observed yields are relatively moderate, approximately 4000 kg/ha. The Maple City output captures many of the years with very high yields, such as 2001 and 2009.  For the northwest region, the Side Green and Bloom outputs are very similar to the output from Maple City, as expected since the Side Green and Bloom input datasets heavily rely on the Maple City climate data. In general, as the Maple City data series was shifted forward in time, more damage occurred as colder temperatures occurred at later dates when the phenology is advanced and therefore more sensitive to cold temperatures. Correspondingly, as the temperatures were shifted backward in time, less damage occurred. The outputs using the Side Green and Bloom data series are very similar to each other, but overall the output using the Side Green data has slightly fewer buds remaining than the output from the Bloom data.  44    Figure 2.7: Model Simulated Buds Remaining and Observed Yields for the Northwest Region Time series of observed yields for the northwest region versus simulated buds remaining. The top left figure is the output from the Maple City climate series, top right shows the output from Side Green and Bloom, bottom left shows the output from Additional Station, and bottom right shows the output from PRISM and NLDAS data.   45  The Additional Station model output has a much shorter series (2005-2015) due to data availability. The output obtained using this series is fairly similar to that obtained using the Maple City series, although the deviations from the observed yields are smaller in the northwest region for some years. For example, in 2010 the Additional Station output represents the observed yield more accurately, while the Maple City output underestimates the remaining buds. The Additional Station output also performs better in 2005 and 2011. However, in 2014 the Additional Station output overestimates the buds remaining. Overall, the output from the Additional Station series is more representative of the observed yields than the output from the Maple City series for the overlapping period of the series.  For the northwest region, the output obtained from the NLDAS input series shows buds remaining of 1.0 for all years except 1986 where buds remaining decreased to 0.98. This model did not accurately represent the variability of the observed yields. The output from the PRISM input has more variability than NLDAS, but buds remaining only ranges from 0.8 to 1.0 in all but three years. In 1986 buds remaining decreases to 0.57, in 1987 to 0.75, and in 2012 buds remaining decreases to approximately 0.3. Overall, the PRISM output does not display as much variability as the observed yields.  The spatially weighted model output performs similarly in the west central region as it did in the northwest (Figure 2.8). The output from the Hart series accurately represents the low yields that occurred in 1986 and 2012, but misses the low yields in 1983 and 2012. In general, the Hart output overestimates the buds remaining, showing high buds remaining when the observed west central yields are low. The Side Green and Bloom output again are similar to each other, but differ from the output obtained from the Hart series because they are based on the adjusted Maple City series. For example in 1983 and 2002, the low observed yields are better  46    Figure 2.8: Model Simulated Buds Remaining and Observed Yields for the West Central Region Time series of observed yields for the west central region versus simulated buds remaining. The top left figure is the output from the Hart climate series, top right shows the output from Side Green and Bloom, bottom left shows the output from Additional Station, and bottom right shows the output from PRISM and NLDAS data 47  represented by the Side Green and Bloom output than the output from the Hart series. Additionally, the Side Green and Bloom outputs have greater variability than the Hart output. As seen in the northwest region, the Side Green output has slightly lower buds remaining values than the Bloom output in general. The Additional Station output is more representative of the observed yields than the Hart output for the west central region. In contrast to the representativeness of the Additional Station output, values for the NLDAS output equal 1.0 for all years. The PRISM output has larger variability than the NLDAS output, and represents the observed yields well in some years such as 2010, but has less variability than the observed yields in the 1990s.  In this analysis, the observed yields are regional while the sour cherry model was developed for a single point location. The representativeness of a single station for the region may be poor due to microclimate influences, and the results may have differed if other locations were used. This may account for the poor performance of the simulated buds remaining using the Hart time series for the west central region.  Overall, the output from additional point data is more representative of the observed yields than that from the single point and gridded climate data variants, indicating the importance of additional data in estimating yields over large areas. Although gridded climate data provide a spatially continuous data surface, NLDAS less suitable in the sour cherry model application, possibly due to the low resolution and the close proximity to Lake Michigan which may have affected the values of the grid cells. The PRISM output performs better in the west central region than the northwest region, but is not as accurate as the Additional Station output in both regions. This is likely a result of the limited number of original input data sites available for the PRISM data series. In addition, while topography is included in the PRISM data development, the 4 km 48  spatial resolution of the data set is generally too coarse to represent the smaller microclimatic variations across the study region. The Additional Station variant performs the best because it contains more basic input information and better captures the spatial variation across the production regions. Therefore, investigators should always carefully consider the input data series used in the development of gridded datasets, especially in applications where microclimate is known to play a role. Finally, it is important to note that these sour cherry simulation results do not account for any other yield-limiting factors such as insects, weeds, or diseases, or weather conditions during pollination. The results should thus be interpreted as patterns or trends under idealized conditions.    49  CHAPTER 3: HISTORICAL AND FUTURE YIELDS 3.1 Background 3.1.1 Susceptibility of Perennial Crops to Climate Change As temperatures warm due to climate change, there is uncertainty how perennial fruit crops will be impacted. In temperate areas, a warming climate may lead to advanced phenological development causing plant tissue to be vulnerable earlier, increasing possibility of frost injury (Pagtera and Arorab, 2013). Other factors such as pollination, disease, and insects also influence fruit production, but in temperate climates springtime freezes are the most limiting factor (Winkler et al., 2013). Additionally, under a warming climate the climatological spatial limits of perennial production may shift as most tree fruit species have a minimum temperature threshold beyond which tree injury occurs (Daly et al., 2012). For example, sour cherry trees can be damaged by temperatures colder than -34°C (Dennis and Howell, 1974). Increasing temperatures may pose a problem for perennial crops in already warm climates like California, as the chilling requirement may not be met in the future (Luedeling et al., 2009). Perennial fruits are especially susceptible to climate change as the capacity for adaptation is limited. It takes 15 to 30 years to develop new perennial crops with selective breeding or genetic engineering, and additional years for the crop to reach production potential (Hatfield et al., 2014).  Assessing the impact of climate change on perennial fruit production is complex and involves consideration of the uncertainty associated with modeling of the climate and its relationship with plant growth and the resulting fruit production. The timing, frequency, and severity of springtime freezes are critical, but changes in springtime freezes are difficult to project. With warmer temperatures phenology may advance earlier, causing increased risk for frost damage, but the timing and frequency of the freezes may also change, causing little change 50  or decreased risk (Fitchett et al., 2014; Winkler et al., 2013). Therefore, it is uncertain whether springtime frost risk for fruit crops will increase or decrease in temperate regions in the future.  3.1.2 Historical Trends in Phenology and Frost Risk Warming temperatures have impacted phenology and frost risk of plants globally. Many studies have reported an advancement in phenology over the past few decades due to the warming temperatures. For example, Menzel et al. (2006) established that phenology has advanced for hundreds of plant varieties, based on an extensive phenology observation network in Europe. Historical trends in the timing of bloom have also been examined, especially for fruit trees. Chmielewski et al. (2004) determined that the start of development and timing of blossom for sweet cherries in Germany advanced by one week in the past 40 years, mainly due to warming in February and April. Similarly, sour cherry bloom in central Europe advanced three days per decade due to increasing April temperatures (Kurlus et al., 2013). A study of apple bloom using four different phenology models concluded that bloom has been advancing in Europe, but more rapidly in the continental areas in Western Europe than the Mediterranean climate regions (Legave et al., 2013). Additionally, phenological observations indicated that flowering across multiple cultivars for apple and pear trees in France and Switzerland advanced due to warming in February and March; however, the changes did not impact fruit production (Guédon and Legave, 2008).  While warming temperatures have caused advances in phenology, the spring frost risk also depends on the timing, frequency, and severity of sub-freezing temperatures. The timing of the last frost compared with a phenological stage is an indicator for the frost risk, as a frost occurrence can be detrimental to tree buds in late phenological stages. Fewer studies have focused on the historical trends in frost risk, and these studies show varying results. Dai et al. 51  (2013) determined that frost risk for woody plants in eastern China decreased over the past 50 years. They concluded that frost risk decreased because phenology advanced by 1.9 days per decade while the date of the last frost advanced at approximately the same rate or at a slightly faster rate of  2 days per decade. Additionally, the number of frost days decreased during the study period.  Similarly, Scheifinger et al. (2002) investigated frost risk at 50 climate stations in central Europe and determined frost damage risk decreased from 1951 to 1997, based on phenological data for 13 plant stages. Some of the plant species included in their analysis were not necessarily late frost sensitive plants, but were used to illustrate the relationship between frost and phenology nevertheless. A frost was defined as a temperature less than -1°C. General warming caused an advancement of phenology and last frost dates. The last frost occurrence advanced by 0.2 days per year and phenology advanced by a range of 0 to 0.2 days per year, with the last phenological stage advancing at a faster rate of -0.28 days per year. Because the last date of frost advanced faster than most phenological stages, they concluded that the risk of late frost damage for earlier phenological stages had decreased.  In contrast, frost damage to woody species in Illinois (U.S.) increased in recent years due to advancements in phenology with warmer March temperatures but little change, or a slight increase, in frequency of April temperatures -1.7°C (Augspurger, 2013). Additionally, Fitchett et al. (2014) found that the frost risk of citrus fruits in Iran had increased historically. Although both the peak flowering dates and the dates of last frost advanced over time, the date of last frost advanced slower than the phenology, increasing the risk of frost damage. They argued that the differences in rates are due to the influence of both maximum and minimum temperatures on the time of flowering, while the occurrence frost events is only a function of minimum temperatures. 52  Warmer minimum temperatures during the study period did not reduce the frost risk, but did lessen the severity of the frost.  In sum, although there is overall agreement on the advancement of phenology due to warming temperatures, trends in frost risk appear to vary by location.  3.1.3 Projected Future Changes in Phenology and Frost Risk As the climate continues to warm, the vulnerability of fruit crops to springtime freeze damage may change. Several studies have projected advancements in phenology under climate change. For example, several statistically downscaled Global Climate Models (GCMs) under differing emissions scenarios projected earlier peach and pear flowering dates in South Korea by the end of the century (Hur and Ahn, 2014). With advanced phenology, freeze risk may change depending if the timing of spring freezes also becomes earlier or if the frequency of freezes decreases. Several studies have addressed this concern.  Ladányi et al. (2010) investigated frost risk for sour cherries in Hungary and concluded risk may decrease by the mid-century. They used the historical mean bloom length as a reference to examine frost risk for 2021-2050 based on the dynamically-downscaled RegCM3.1 Regional Climate Model (RCM) scenario. For the future period, no frost days were projected to occur during bloom or 10 days prior to the beginning of bloom, due to significant increases in the mean minimum temperature. Additionally, an increased number of days with precipitation during bloom was projected, which indicates a worsening of conditions for pollination.  Decreasing frost damage risk is also projected for viticulture in Luxembourg. Molitor et al. (2014) used six dynamically-downscaled RCM simulations that were further empirically downscaled and bias corrected using quantile mapping. To determine future frost risk, the date of budburst was estimated using a chilling model and a phenological forcing model with a 53  photoperiod variable. A frost occurrence within 60 days after budburst was considered to be damaging. The median date of budburst for the end of the century advanced by 11 days, and the date of last frost advanced by 28 days relative to the historical period. Therefore, the length between the date of last frost and day of budburst increased, indicating less frost damage in the future.  Similarly, Eccel et al. (2009) examined the risk of spring frost to apples in Italy for historical and future periods by modeling break of dormancy and phenology using growing degree days (GDDs). For the historical period, the flowering date at one station advanced, on average, by 0.88 days per year and at the other station by 0.67 days per year. In spite of the earlier phenology, a decrease in frost risk for apples was seen in recent decades. Extending the analysis to a future period, they used linear regression to statistically downscale simulations from the HadCM3 GCM simulation under multiple emissions scenarios to the location of the two stations. Flowering from 1991 to 2056 was projected to advance, ranging from 0.04 to 0.13 days per year. The range of uncertainty was larger for future frost risk, with 11 scenarios suggesting decreased risk, 6 scenarios showing increased risk, and 7 scenarios showing no change. In general, the frequency of frost episodes was larger for the higher emissions scenario, Special Report on Emissions Scenarios (SRES) A2, and smaller for the lower emissions scenario, SRES B2. Overall, the authors concluded that frost risk is likely to stay the same, or slightly decrease, in the future. Hoffmann and Rath (2013) used 7 phenology models and simulations from 13 RCMs under emissions scenario SRES A1B that were bias corrected using quantile mapping to investigate future frost risk for apples in Lower Saxony, Germany. On average the last spring freeze advanced 10 days by 2035 and 27 days by 2084. Generally, phenology advanced slower, 54  increasing the time between bloom and the last spring freeze, indicating decreased frost risk in the future. However, one scenario, which projected the largest advancement in bloom, did project an increase in frost risk.  Chmielewski et al. (2010) also looked at frost risk for apple trees in 11 fruit growing regions in Germany, but, in contrast to Hoffmann and Rath (2013), found increasing frost risk in the future. The ECHAM5/OM GCM under SRES A1B was statistically downscaled using WETTREG (weather situation-based regionalization method) for the future period of 2011 to 2100. Five phenology models were used to simulate apple blossom, and the frequency of late frost events up to 10 days after the beginning of apple blossom were examined. A temperature threshold of -2°C was used to represent light frosts which would harm the flower buds by 10%, and -4°C for a damage of 90%. Blossom advanced about 15 days by 2100, and the mean probability for late frost after the beginning of blossom increased by 8%. Generally, the frequency of light frosts increased, and the authors concluded that the increase in damage risk will most likely be moderate in the future. Mosedale et al. (2015) examined frost risk for grapevines in southwest England and also supported an increase in future frost risk. An ensemble of climate projections from RCMs driven by the HAdCM3 GCM for low, medium, and high emissions scenarios were stochastically downscaled with a weather generator, producing 1000 weather sequences. Phenology was modelled using GDDs with various fixed starting dates and also with a simulated break of dormancy using winter chilling requirement. They evaluated frost risk as the probability of frosts after budbreak, although results were similar using measures capturing the severity or frequency of frost events. By 2080 the scenarios projected an advancement in the modelled date of budbreak by 45 days, an advancement in the date of last frost by 19 days, and therefore an 55  increased frost risk. However, projected frost risk was highly sensitive to the choice of phenology model and to the simulated start of budbreak. Additionally, adverse weather such as precipitation or low temperatures during flowering declined during the future period, indicating improved flowering conditions.  Future change in frost risk may also depend on the location as shown by Kaukoranta et al. (2010) for 9 apple cultivars at 14 locations in Finland. A simulation from the RCA3 RCM driven by the ECHAM4/OPYC3 GCM under SRES A2 was downscaled using the delta method for the near future, 2011-2040. They simulated the timing of flowering using accumulated GDDs from January 1, and examined the frequency of temperatures less than -2°C after flowering. While frost risk was generally expected to not change, slight decreases in the west and increases in the south were projected with the early cultivars being particularly vulnerable to future frost damage.  Rochette et al. (2004) also concluded that frost risk varied by location. They used simulations under three emissions scenarios from a single GCM (CGCMI) for two future periods, 2010-2039 and 2040-2069. The projections were downscaled using the delta method for 69 locations in Canada. Apple tree phenology was modelled using accumulated GDDs starting on January 1. Earlier dates of last frost and phenology were projected for all regions but the rates of advancement differed. In eastern Canada, the frost risk was unchanged, but the projections suggest small increases in frosts risk in southern Quebec, large increases in southern Ontario, and large decreases in frost risk in the continental north. The large increases in frost risk in southern Ontario were due to larger increases in temperature during winter than at other locations as well as the already mild winter climates in southern Ontario compared to other regions, causing a faster rate of phenology advancement.  56  Winkler et al. (2012b) used a large ensemble of climate projections to examine freeze risk for sour cherry trees in Michigan. The projections were developed using multiple downscaling approaches from four GCMs driven by two emissions scenarios (SRES A2 and B2). Daily GDD accumulation beginning on January 1 was used to estimate bud development, and the last occurrence of freezing temperature was used to estimate freeze risk. Phenology was projected to advance in the future, and the date of the last spring freeze was projected to advance as well, but with a larger uncertainty range. A large number of scenarios projected little or no change in freeze risk in the future, whereas a similar number of scenarios suggested either an increase or decrease. Therefore, they concluded that the uncertainty surrounding spring freeze risk is large.  The studies summarized above demonstrate the considerable uncertainty about future frost risk and the need for more research to understand climate change impacts on perennial crop production. The differing results of the previous studies also highlight the importance of using an ensemble, or multiple climate projections, as the projected changes can vary depending on choice of climate projection. While some studies estimated future changes in frost risk, none were able to transform the projections of future frost risk into fruit production loss. Additionally, spatial coverage of previous analyses is uneven. In particular, few studies were conducted for the temperate climate regions in the U.S. 3.1.4 Climate Projections Many climate change impact studies use an ensemble of climate projections. The choice of climate model, emissions scenario, and downscaling method can lead to varying future projections (Hanssen-Bauer et al., 2005). Ensembles are generally used because there is no one the ensemble will provide an uncertainty range in the climate projections 57  (CCSP, 2008; Winkler et al., 2011b; Winkler et al., 2012a). However, different GCMs or downscaling methods may be similarly biased, so the ensemble may not capture the entirety of the uncertainty (Hanssen-Bauer et al., 2005).  Climate projections generally need to be downscaled for applications as GCMs typically have coarse spatial resolutions, ranging from 100 to 300 km (Winkler et al., 2011a). Downscaling methods can be described as empirical or dynamical. Empirical downscaling derives regional climate information from a statistical model that relates large scale climate variables, or predictors, to more local variables, or predictands (Wilby et al., 2004). Climate model output is fed into the statistical model to estimate the local climate characteristics, with the major assumption that the statistical relationships hold true under the future climate (Wilby et al., 2004). The most basic and commonly used form of statistical downscaling is the change factor, or delta, method. In delta downscaling, the difference between the GCM projected values for a future and a control period is applied to a historical time series. However, when using this method the range and variability of the future climate do not change relative to the historical climate (Wilby et al., 2004; Winkler et al., 2012a). Dynamical downscaling is when numerical models are used to simulate a finer resolution climate projection, like nesting a RCM within a GCM. RCMs have finer spatial resolutions, approximately 25-50 km, and are typically used when regional or local influences, such as topography, play important roles in the regional climate (Winkler et al., 2011a).  Climate projections often need to be adjusted for biases, and it can be difficult to separate bias correction from downscaling (Winkler, 2016). The choice of correction is important as climate projections are sensitive to bias correction methods (Watanabe et al., 2012). Quantile mapping (QM) is a frequently used downscaling and bias correction method. QM corrects for 58  distribution bias of the climate projections via comparisons of the observed and modeled cumulative probability distributions. Multiple linear regression (MLR) is another widely used downscaling and debiasing approach. A linear function between predictors and a predictand is developed and solved so that error is minimized. Multiple linear contour regression (MLCR) developed by Abraham et al. (2014) is a variant of MLR that minimizes the prediction error in the climate projections and removes biases in the predicted distribution through an iterative process based on the solution of MLR. This method showed improvement over standard regression methods for maximum and minimum temperatures and precipitation for climate stations in Michigan (Abraham et al., 2014). All empirical downscaling and debiasing methods have strengths and weaknesses. For example, whereas QM can reduce biases of daily temperatures and precipitation (Themeßl et al., 2012), the prediction error is higher compared to linear regression method (Abraham et al., 2014). On the other hand, MLR does not correct the distribution bias, and is therefore not recommended for precipitation error correction (Themeßl et al., 2011). 3.1.5 Climate Change Impacts on Perennial Crops in Michigan Over 70 percent of U.S. sour cherry production occurs in Michigan. The Great Lakes are critical for sour cherry production in Michigan, modifying the local climate and reducing the risk of winter kill and spring frost. Assessing historical and future frost risk for this region is challenging, however. In the Midwest, between 1900 and 2010 the average temperature increased more than 0.83°C (1.5°F), with the highest rate of increase since 1980 (Kunkel et al., 2013). Generally, the greater proportion of warming for 1895-2010 occurred during the winter and spring seasons, with much of the warming associated with increased minimum temperatures (Andresen et al., 2012). While winters have warmed and the last spring freeze is earlier, the 59  interannual variability of extreme cold outbreaks remains large (Kunkel et al., 2013). Additionally, precipitation has increased during the past century in the Midwest, mostly due to more intense rainfall events (Pryor et al., 2014), which may have influenced pollination conditions. Furthermore, the spatial extent of plant hardiness zones, which are based on annual extreme minimum temperatures, has shifted northward one half zone in Michigan (Daly et al., 201).  Temperatures are projected to increase in the Michigan sour cherry production region. In the Midwest, temperatures are projected to increase by 2.1-2.7°C by the mid-century and by 3.1-4.7°C by the end of the century depending on emissions scenario (Pryor et al., 2013). Both GCMs and RCMs project that the annual number of days below freezing as well as extreme freezes will decrease in Midwest (Kunkel et al., 2013). Correspondingly, cold air outbreaks are expected to decrease in the future across most of the U.S. (Vavrus et al., 2006). Moreover, the date of last spring frost is projected to occur 2-5 weeks earlier by the end of the century in the Midwest (Wuebbles and Hayhoe, 2003). Precipitation is projected to increase during winter and spring on average in the Great Lakes Region (Hayhoe et al., 2010). However, models generally have larger uncertainty with precipitation than temperature (Pryor et al., 2014). Additionally, changes in climate variability may significantly impact agriculture, but uncertainty surrounds the future changes in the variability of rainfall and temperature (Thornton et al., 2014). Climate projections for Michigan are complicated, however, by the poor representation of the Great Lakes in climate models. The coarse resolution of GCMs does not account for the Great Lakes, and many RCMs do not include lake models, instead estimating the lake temperatures as an average of the Atlantic and Pacific Ocean temperatures near the coast (Winkler et al., 2012a). Neither the RCMs nor GCMs consistently control the temperature variability in the Great Lakes 60  Region, reflecting the inconsistent treatment of the Great Lakes in climate projections (Mearns et al., 2013). 3.1.6 Study Objectives  The objectives of this study are to examine how sour cherry yield has changed historically in Michigan, how climate change may impact yield in the future, and to evaluate how the choice of climate projection can influence the interpretation of future changes. This study uses an ensemble of 88 individual scenarios obtained from various climate models, emissions scenarios, and downscaling methods available for three locations in the fruit growing regions of the western Lower Peninsula of Michigan. The climate projections serve as input to a sour cherry yield model which utilizes relationships between climate, phenology, freeze damage, and pollination conditions to estimate yield.  61  3.2 Data and Methods 3.2.1 Sour Cherry Yield Model A sour cherry yield model developed by Black et al. (in preparation) was used in this study to examine both historical and future trends in yield. This hybrid empirical/process-based model simulates eight stages of phenology based on the work of Zavalloni et al. (2006) using summed base 4ºC GDD thresholds. The GDDs accumulate beginning on January 1, under the assumption that the chilling requirement will be met by that date under projected future conditions. The GDD thresholds used in Chapter 2 were adjusted accordingly by 10 GDDs for the earlier break of dormancy date. Freeze damage was estimated using critical temperature thresholds for each phenological stage based on a chamber study by Dennis and Howell (1974). The cumulative seasonal damage was subtracted from the buds remaining variable which is a ratio from 0 to 1, indicating survival of the sour cherry buds. The model also accounted for the effects of poor pollination on yield. Poor pollination conditions, when bee activity is not favorable, was represented as a summed number of days when precipitation occurred or the average temperature was less than 10°C during phenological stage 8 (full bloom). A yield algorithm used both the cumulative freeze damage and frequency of poor pollination days to estimate sour cherry yield per year. The theoretical potential yield was 14308 kg/ha which would occur during a year without any damaging freezes or poor pollination days. This model was validated against observed yields aggregated for the northwest region of Michigan from the United States Department of Agriculture (USDA) National Agricultural Statistics Service (NASS). Simulated yield at Maple City accounted for 31% of the variation of the observed yield and was statistically significant at the 0.001 alpha level. The sour cherry yield model was used to assess historical and future temporal trends in yield, buds remaining, poor 62  pollination days, phenology, and damaging freeze events. The model output from a historical reference period 1980-2000 was used as a baseline to evaluate future changes in sour cherry yield. 3.2.2 Historical Climate Data Historical climatic data series of daily maximum and minimum temperature and precipitation were obtained from the National Weather Service (NWS) Cooperative Observer Program (COOP) for one reference site in each of the three major sour cherry production regions, namely, Maple City (northwest), Eau Claire (southwest), and Hart (west-central), Michigan. These observations were used in the sour cherry model to examine historical temporal trends for 1960-2015. The data series at Hart extended to 1894, and these additional years were also analyzed. Additionally, average annual temperatures for northwest Michigan, the third climate division, for 1960-2015 from National Centers for Environmental Information (NCEI) Climate at a Glance were used to examine historical temporal trends in temperature. Because it is common to move climate stations and/or update equipment, the climate observations for the historical reference period (1980-2000) were tested for homogeneity using the methodology by Wijngaard et al. (2003). The seasonal and annual average maximum and minimum temperatures and precipitation were evaluated with the standard normal homogeneity test for a single break (Alexandersson, 1986), the Buishand range test (Buishand, 1982), the Pettitt test (Pettitt, 1979), and the Von Neumann ratio test (Von Neumann, 1941). All three climate series passed the majority of the tests. 3.2.3 Sources of Future Climate Projections To examine future trends in sour cherry yields, multiple climate projections developed by CLIMARK, a National Science Foundation (NSF) funded project focused on an integrated 63  framework for climate change impact assessments for international market systems with long-term investments (http://cherry.geo.msu.edu), were used in the sour cherry yield model. These projections included output from 16 GCMs from the 5th phase of the Coupled Model Intercomparison Project (CMIP5) (Table 3.1). These models were chosen as simulations for three representative concentration pathways (RCPs). The resolution of the CMIP5 models ranges from 0.5° to 4° for the atmosphere component and from 0.2° to 2° for the oceanic components, depending on climate model (Taylor et al., 2012). This study used the long-term (century time scale) CMIP5 experiments which were driven by RCPs 4.5, 6.0, and 8.5 which show a range of radiative forcing and concentration outcomes based on socioeconomics and mitigation strategies (Moss et al., 2010). The historical CMIP5 simulations were used as a control period.  Additionally, dynamically-downscaled simulations available from the North American Regional Climate Change Assessment Program (NARCCAP) were used (Table 3.2). A total of.8 simulations obtained from a combination of 4 RCMs and 4 GCMs were available (Table 3.3). The NARCCAP simulations have a 50 km spatial resolution (Mearns et al., 2007, 2009). The simulations were forced with the SRES A2 emissions scenario for the 21st century. The A2 emissions scenario represents a very heterogeneous world based on socioeconomic factors such as increasing population, relatively slow technological change, and regional economic growth (Nakicenvoic et al., 2000). NARCCAP simulations were available for a control period (1978-1998) and for a mid-century period (2040-2060).  The Great Lakes were crudely represented, if at all, in the CMIP5 models. All of the RCM included realistic outlines of the Great Lakes, but the methods used to estimate lake temperatures differed. The Canadian Regional Climate Model (CRCM) contained a thermodynamic mixed layer lake model for the Great Lakes Region (Govette et al., 2000). In   64  Table 3.1: Global Climate Models used in the CMIP5 Projections The CMIP5 GCMs used in the construction of climate projection ensemble. Modeling Center Model Name Institution FIO FIO-ESM The First Institute of Oceanography, SOA, China IPSL IPSL-CM5A-LR  IPSL-CM5A-MR Institut Pierre-Simon Laplace, France MIROC MIROC5  (Model for Interdisciplinary Research on Climate) Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology MIROC MIROC-ESM  MIROC-ESM-CHEM Japan Agency for Marine-Earth Science and Technology, Atmosphere and Ocean Research Institute (The University of Tokyo), and National Institute for Environmental Studies MOHC  (additional realizations by INPE) HadGEM2-ES Met Office Hadley Centre (additional HadGEM2-ES realizations contributed by Instituto Nacional de Pesquisas Espaciais) MRI MRI-CGCM3 Meteorological Research Institute, Japan NASA GISS GISS-E2-H  GISS-E2-R NASA Goddard Institute for Space Studies, USA NCAR CCSM4 National Center for Atmospheric Research, USA NIMR/KMA HadGEM2-AO National Institute of Meteorological Research/Korea Meteorological Administration NOAA GFDL GFDL-CM3  GFDL-ESM2G  GFDL-ESM2M NOAA, Geophysical Fluid Dynamics Laboratory, USA NSF-DOENCAR CESM1 (CAM5) National Science Foundation, Department of Energy, National Center for Atmospheric Research, USA  65    Table 3.3: Available NARCCAP Simulations Available RCM and GCM combinations for NARCCAP simulations.  RCM GCM-driven NCEP-driven   GFDL CGCM3 CCSM HADCM3   RCM3 X X     X CRCM   X X   X WRFG   X X   X HRM3 X     X X   Table 3.2: Climate Models used in the NARCCAP Projections  The regional and global climate models used to produce the NARCCAP dynamically-downscaled climate projections. Modeling Type Model Name Institution Regional Climate Model CRCM  Canadian Regional Climate Model Canadian Centre for Climate Modelling and Analysis Regional Climate Model HRM3  Hadley Regional Model 3 Met Office Hadley Centre, UK Regional Climate Model RCM3  Regional Climate Model Version 3  (also referred to as RegCM3) International Center for Theoretical Physics, Italy Regional Climate Model WRFG  Weather Analysis and Forecast  Model: Grell Convection Scheme National Center for Atmospheric Research, USA Global Climate Model CCSM  Community Climate System Model National Center for Atmospheric Research, USA Global Climate Model CGCM3 Coupled Global Climate Model 3 Canadian Centre for Climate Modelling and Analysis Global Climate Model GFDL Coupled Model CM2.1 NOAA, Geophysical Fluid Dynamics Laboratory, USA Global Climate Model HADCM3 Hadley Center Coupled Model Version 3 Met Office Hadley Centre, UK   66  contrast, the Weather Analysis and Forecast Model (WRFG) simulation driven by the Community Climate System Model (CCSM) assigned the lakes temperatures based on the GCM surface temperatures; however, the CCSM model did contain a thermodynamic lake model that simulated surface temperatures for larger lakes (UCAR, 2007). The Regional Climate Model Version 3 (RCM3) driven by the Coupled Model CM2.1 (GFDL) calculated the lake temperatures by interpolating the ocean sea surface temperatures from the GCM to the same latitudes over the Great Lakes (UCAR, 2007). For the RCM3 model driven by the Coupled Global Climate Model 3 (CGCM3), the lake temperatures were assigned by extrapolating the GCM temperatures downward, which was generally only 2 data points for the Great Lakes (UCAR, 2007). For the Hadley Regional Climate Model (HRM3) driven by GFDL, sea surface temperature ancillary data from the HadCM3Q0 GCM were applied to the GFDL data to simulate temperatures over the Great Lakes, which caused issues with the resulting projection having cooler future lake temperatures in the summer than the current conditions (UCAR, 2007). Time series of maximum and minimum temperature and precipitation at a daily time step were extracted from the CMIP5 archive. Daily maximum and minimum temperature files were also available for the NARCCAP simulations. The NARCCAP precipitation outputs had a three-hourly temporal resolution, which were aggregated to obtain a daily precipitation total. 3.2.4 Downscaling and Debiasing Methods Twenty-one year control and future periods were employed. The definitions of the control and future periods were contingent up on the length of the time series available for the CMIP5 and NARCCAP simulations. For CMIP5, the control period was defined as 1980-2000, whereas for NARCCAP the control period was 1978-1998. The future period was defined as 2040-2060 for both simulation sources.  67  The climate projections were downscaled to the three Michigan locations, Eau Claire, Hart, and Maple City, using various methods (Table 3.4). For the CMIP5 scenarios, the delta, or change factor, method was employed. For the NARCCAP scenarios, the downscaling methods of delta, quantile mapping, multiple linear regression, and multiple linear contour regression were used. Additionally, time series were taken directly from the nearest grid point to a station, and these data will be referred to as the raw NARCCAP projections.  Table 3.4: Downscaling and Debiasing Methods Downscaling and debiasing methods and the baseline data used to calculate future changes. Downscaling Method Abbreviation Description Baseline Delta approach Delta Monthly change factors are computed between future and control climate simulations and are applied directly to observations Observations  (1980-2000 for CMIP5,  1978-1998 for NARCCAP) Raw, or pure, downscaling Raw Extracts data from the nearest grid point to a station location Simulated control  (1978-1998) Quantile mapping QM Adjusts the distribution of the climate simulations to match the distribution of the observations Simulated control  (1978-1998) Multiple linear regression MLR Minimizes the error of the climate simulations by solving a linear regression between observations and control simulations, and applying the corrections to the future simulations Simulated control  (1978-1998) Multiple linear contour regression MLCR Simultaneously minimizes the prediction error and adjusts the distribution of the climate simulations, and the corrections are applied to the future simulations Simulated control  (1978-1998)   The delta downscaling method is simple statistical downscaling where monthly change factors between the future and control climate simulations were calculated and applied to temperature and precipitation observations for a historical period (1980-2000). The change 68  factors were a simple difference for maximum and minimum temperatures which were added to the temperature time series for the historical period to obtain daily time series for the future period. The change factors were defined as a ratio for precipitation, and the historical precipitation time series were multiplied by the change factor to obtain future precipitation series. Quantile mapping (QM) was also used for the NARCCAP RCM/GCM scenarios (Liszewska et al., 2012). QM is a quantile-based empirical error correction method (Themeßl et al., 2012). This method corrected the shape of the control run distribution of the simulated daily time series of the climate variables, maximum and minimum temperature and precipitation, based on the distribution of the observations, and the corrections were then applied to the future simulations. The downscaling method of multiple linear regression (MLR) was also used. MLR is a commonly used downscaling method that uses a linear regression to model the relationship between explanatory variables and response variables, and minimize the error of the downscaled values (Themeßl et al., 2011). In this study, the functions were developed seasonally using daily maximum and minimum temperature and precipitation observations as the response variables, and a set of surface and upper-air climate variables as the explanatory variables. The functions were developed using time series of the explanatory variables from the RCM simulations driven by a global reanalysis and observations of the responsible variables. The regression functions were then applied to the future climate simulations to generate projections for 2040-2060. Moreover, multiple linear contour regression (MLCR) developed by Abraham et al. (2014) was used as it minimizes the prediction error and removed biases in the predicted distribution. This method used a linear regression and a squared error loss function to model the relationship between the observations of the response variables and reanalysis-driven simulations of the 69  explanatory variables, while simultaneously correcting the distribution of the control simulations with the observed distribution. The corrections were then applied to the future projections.  3.2.5 Application of the Climate Projections  The future and control climate projections of daily maximum and minimum temperature and precipitation served as input to the sour cherry yield model. The average model output for 2040-2060 was compared to: 1) for the delta projections, the model output when the input climate series were the observed maximum and minimum temperature and precipitation series l output when the input climate series were the RCM-simulated maximum and minimum temperature and precipitation time series for the control period, and 3) for all other climate projection times, the model output when the input climate series were the downscaled and debiased maximum and minimum temperature and precipitation time series for the control period. This analysis was performed for the variables yield, buds remaining, poor pollination days, date of stage 2, and date of stage 8. Additionally, the frequency of damaging freeze events over the 21-year period and the average severity of damage per damaging freeze were examined. Percent changes for yield, buds remaining, and average severity were calculated by subtracting the average baseline values from the average future value, then dividing by the average baseline value, and multiplying by 100. Further analyses were completed to examine why yield increased or decreased, using seasonal changes in climatic variables.  70  3.3 Results 3.3.1 Historical Temporal Trends Linear trends in simulated yield for the period 1960-2015 are statistically insignificant for the three stations, Eau Claire, Hart, and Maple City (Table 3.5), but inspection of the time series of simulated yield (Figure 3.1) suggests non-linearity in the magnitude and direction of the trend at Hart, although not necessarily at the other two stations. Simulated yields increase at Hart from 1960 until the 1990s, but decrease in the last 2 decades. Similar complex trends are seen for the extended period (1894-1959) at Hart, with a decrease in simulated yield until the 1940s followed by an increase until 1959. In general, the simulated historical yields are smaller for Maple City than Eau Claire and Hart. In general, the simulated values of buds remaining are higher at Eau Claire and Hart than at Maple City (Figure 3.1). Although the simulated temporal trends in buds remaining for the historical period and the differences among the three stations visually correspond with those for simulated yield, the linear trends, like those for yield, are statistically insignificant at all three stations (Table 3.5). Temporal trends in the estimated number of poor pollination days are also insignificant, although some nonlinearity is evident. For example, the time series for Hart shows a general decrease in the number of poor pollination days from 1960 to the 1980s, followed by an increase in number of poor pollination days to the present, whereas the simulations for the additional years of record at Hart suggest little change or an increase in poor pollination days. At Maple City, the estimated number of poor pollination days decreases from 1960 to 2000, followed by an increase until the present. In contrast, a general, but insignificant, increase in the estimated number of poor pollination days is evident at Eau Claire.  71     Table 3.5: Linear Trends in Model Output at Eau Claire, Hart, and Maple City, 1960-2015 Historical linear trends by variable for each station, Eau Claire, Hart, and Maple City for the period 1960-2015. The linear trends are characterized by slope, p-value, and R-squared.    Station Slope P-value R2 Yield Eau Claire -26.63 0.326 0.02 Hart -5.53 0.836 0.00 Maple City 34.96 0.256 0.02 Buds Remaining Eau Claire 0.00 0.739 0.00 Hart 0.00 0.981 0.00 Maple City 0.00 0.100 0.05 Poor Pollination Days Eau Claire 0.02 0.246 0.02 Hart 0.01 0.804 0.00 Maple City 0.00 0.865 0.00 Stage 2 Eau Claire -0.22 0.012** 0.11 Hart -0.06 0.447 0.01 Maple City -0.11 0.188 0.03 Stage 8 Eau Claire -0.19 0.009*** 0.12 Hart -0.01 0.839 0.00 Maple City -0.09 0.205 0.03 Frequency of Damage Events Eau Claire 0.03 0.060* 0.06 Hart 0.01 0.666 0.00 Maple City -0.02 0.477 0.01 Average Severity of Damage Eau Claire 0.00 0.854 0.00 Hart 0.00 0.876 0.00 Maple City 0.00 0.353 0.02 Significance: * for 0.10, ** for 0.05, *** for 0.01  72   Figure 3.1: Time Series of Simulated Yield, Buds Remaining, and Poor Pollination Days at Eau Claire, Hart, and Maple City, 1960-2015 Yield is given in the top row, buds remaining in the middle row, and poor pollination days in the bottom row. The left column is for Eau Claire, middle column for Hart, and the right column for Maple City. The black line is a 9-year moving average.  73  The estimated dates of stage 2 occurrence are similar for Hart and Maple City, but are much earlier at Eau Claire due to its southern location and consequently warmer temperatures (Figure 3.2). The estimated date of stage 2 (side green) advances at Eau Claire and Maple City throughout the period of record at the rates of 2 days per decade and 1 day per decade, respectively, although the linear trend is significant only at Eau Claire (at the 0.05 alpha level, Table 3.5). At Hart, the linear trend is only 0.6 days per decade and insignificant. Similarly, the linear trend in the estimated advancement of stage 8 (full bloom) is only significant at Eau Claire, with a trend of 1.9 days per decade for 1960-2015.   The simulated damaging freeze events are a function in the model of the timing of phenology and the occurrence of freezing temperatures. In general, Maple City experiences a greater number of damaging freeze events than Eau Claire and Hart (Figure 3.3). The linear trend in the frequency of damaging freeze events is weakly significant (p=0.10) at Eau Claire, and insignificant at the other two stations. Visual inspection of the time series points to limitations in the interpretation of a linear trend, however, especially for Hart and Maple City. At Hart, the number of events decreases until 1990, followed by an increase until the present. Similarly, the number of damaging freeze events initially decreases at Maple City (until the early 2000s), but increases in the past decade. An exceptionally large number of simulated damaging freezes are seen in 2010 and 2012 at both Hart and Maple City.  The severity of damaging freezes is measured in terms of loss of buds, with an index ranging from zero (no damage) to one (all buds destroyed). Many years at Eau Claire and Hart did not have any simulated damage events, which complicates the calculation of the linear trend. Although the linear trends are insignificant at all three stations, the time series shown in Figure 3.3 suggest some non-linear variations. At Eau Claire, the simulated severity of damaging74    Figure 3.2: Time Series of Simulated Dates of Stage 2 and Stage 9 at Eau Claire, Hart, and Maple City, 1960-2015 The date of stage 2 is given in the top row and date of stage 8 in the bottom row. The left column is for Eau Claire, middle column for Hart, and the right column for Maple City. The black line is a 9-year moving average. 75   Figure 3.3: Time Series of Simulated Damaging Freeze Events and Average Severity of Damage at Eau Claire, Hart, and Maple City, 1960-2015 The frequency of damaging freeze events is given in the top row and average severity of damage per damaging freeze event in the bottom row. The left column is for Eau Claire, middle column for Hart, and the right column for Maple City. The black line is a 9-year moving average. 76  freezes generally increases until the 1980s and then decreases, whereas at Hart the simulated severity decreases until 1990 after which it increases. At Maple City there is a general, but insignificant, trend of decreasing severity of damaging freeze events over time.  Lastly, the model outputs are examined for the 1980-2000 period, which is the historical period for comparison of the delta (change factor) scenarios described above. For this period, Hart has the largest simulated yield, followed by Eau Claire, and then Maple City (Table 3.6). The spatial variability for buds remaining corresponds with that for simulated yield. For the estimated number of poor pollination days, Eau Claire has the largest value while Hart and Maple City share a slightly lower value. Hart and Maple City have similar average dates of phenological stages 2 and 8, whereas Eau Claire has an earlier average by 2-2.5 weeks. Maple City has the largest number of damaging freeze events during the 21-year period. Eau Claire has approximately half as many events as Maple City, and at Hart the frequency of damaging freezes is small. Contrastingly, the average severity of damage is fairly similar between stations.   Table 3.6: Historical Averages of Yield Model Output, 1980-2000 Historical averages of the yield model output for the 21-year period, 1980-2000, per station. The frequency of damaging freezes is a summed total over the 21-year period. Station Yield (kg/ha) Buds Remaining Poor Pollination Days Date of Stage 2 Date of Stage 8 Frequency of Damaging Freezes Average Damage Severity Eau Claire 8640.3 0.77 3.7 97.3 118.2 29 0.23 Hart 9718.6 0.87 3.1 113.3 133.5 13 0.20 Maple City 7709.7 0.62 3.1 114.3 133.8 50 0.19   77  3.3.2 Projected Future Changes in Climatic Variables  Projected future changes in maximum temperature, minimum temperature, and precipitation by 2040-2060 were examined by simulation source and downscaling method and, if applicable, by RCP. For temperature, the focus was on annual and seasonal means of maximum and minimum temperature. For precipitation, the annual and seasonal frequencies of wet days and average amount of precipitation per wet day were analyzed.   Ensemble (i.e., multi-model) means of projected future changes in the annual values of the climate variables at Eau Claire, Hart, and Maple City are shown in Table 3.7. The ensemble means of the projected changes in the seasonal means and frequencies are generally similar to the annual values, with the exception of slightly larger projected increases in average amount of precipitation per wet day during spring (not shown). Because of these small differences, the discussion below focuses on the differences in the annual means and frequencies.  At all stations, annual mean maximum and minimum temperatures are projected to increase (Figures 3.4-3.5). However, the magnitude of the ensemble means differs by source of climate projection (Table 3.7). (Because of the similarities in the projected changes between stations, figures are only shown for Eau Claire.) As expected given the larger number of simulations, the uncertainty envelope is larger for the downscaled CMIP5 models, with projected changes ranging from 1 to 4.6°C at Eau Claire. The HadGEM2-ES model and the models developed from the MIROC modeling center project larger increases in maximum temperature compared to the other CMIP5 models. The ensemble means are 2.5°C, 2.0°C, and 3.1°C for RCPs 4.5, 6.0, and 8.5 respectively at Eau Claire. The reason for the higher projected change for the RCP 4.5 simulations compared to the RCP 6.0 simulations is that the radiative forcing is larger for RCP 4.5 until approximately the middle of the 21st century, after which the radiative 78   Table 3.7: Multi-Model Mean Changes by the Mid-Century in Temperature and Precipitation Multi-model mean changes by the mid-century in annual mean maximum and minimum temperature, precipitation per wet day and the frequency of wet days for each downscaling method. Station  Downscaling Method Maximum Temperature (°C) Minimum Temperature (°C) Precipitation Amount per Wet Day (mm) Frequency of Wet Days Eau Claire CMIP5 Delta RCP 4.5 2.5 2.5 0.3 0.0 CMIP5 Delta RCP 6.0 2.0 2.0 0.3 0.0 CMIP5 Delta RCP 8.5 3.1 3.1 0.4 0.0 NARCCAP Raw 2.4 2.4 0.3 -8.8 NARCCAP Delta 2.4 2.4 0.2 0.0 NARCCAP QM 2.3 2.3 0.5 -2.1 NARCCAP MLR 1.8 1.7 0.1 2.6 NARCCAP MLCR 2.1 2.0 0.0 -5.9 Hart CMIP5 Delta RCP 4.5 2.5 2.6 0.4 0.0 CMIP5 Delta RCP 6.0 2.0 2.1 0.3 0.0 CMIP5 Delta RCP 8.5 3.1 3.1 0.7 0.0 NARCCAP Raw 2.4 2.4 0.3 -7.4 NARCCAP Delta 2.4 2.4 0.5 0.0 NARCCAP QM 2.4 2.5 0.7 -0.2 NARCCAP MLR 1.7 1.6 0.1 1.3 NARCCAP MLCR 2.0 2.0 0.1 7.6 Maple City CMIP5 Delta RCP 4.5 2.6 2.7 0.5 0.0 CMIP5 Delta RCP 6.0 2.1 2.3 0.5 0.0 CMIP5 Delta RCP 8.5 3.3 3.4 0.7 0.0 NARCCAP Raw 2.3 2.4 0.3 -7.4 NARCCAP Delta 2.3 2.4 0.4 0.0 NARCCAP QM 2.2 2.5 0.5 0.6 NARCCAP MLR 1.7 1.6 0.0 1.4 NARCCAP MLCR 2.0 2.0 0.1 4.0     79     Figure 3.4: Projected Change by 2040-2060 in Average Annual Maximum Temperature (°C) at Eau Claire by Climate Projection 80      Figure 3.5: Projected Change by 2040-2060 in Average Annual Minimum Temperature (°C) at Eau Claire by Climate Projection 81  forcing is higher for RCP 6.0 (Moss et al., 2010).  Direct comparison of the climate projections from CMIP5 and NARCCAP is difficult because of the different driving greenhouse gas emissions. Very generally, the radiative forcing of the SRES A2 emission scenarios falls between that of RCP 6.0 and 8.5 (Walsh et al., 2014). The ensemble means of annual maximum and minimum temperature for the NARCCAP-derived projections are smaller than those for the CMIP5 RCP 8.0 simulations and more similar to those for the CMIP5 RCP 4.5 and 6.0 simulations (Table 3.7). Differences by downscaling method are evident among the NARCCAP-derived projections, with increases in the ensemble means for mean maximum (minimum) temperature ranging from 1.7°C (1.6°) to 2.4°C (2.5°C) across various downscaling methods. Additionally, the uncertainty envelope for the NARCCAP models, when summed across all the downscaling methods, is smaller than the downscaled CMIP5 models, ranging from 1.2 to 3.2°C with, in general, the WRFG_CGCM3 simulation projecting the smallest increase in maximum temperature and the HRM3_GFDL simulation projecting the largest increase.  Considerable variations in the projected changes in the annual average precipitation amount per wet day are evident (Figure 3.6). The majority of the downscaled CMIP5 simulations project an increase in precipitation amount per wet day. The projections range from -0.3 to 1.8 mm per day across the three RCPs, with the climate model GFDL-CM3 projecting the largest increases. Similarly, the majority of the NARCCAP models project increases in precipitation, with the exception of the projections developed with the MLR and MLCR downscaling methods. For these two methods, an equal number of the projections suggest small decreases in    82       Figure 3.6: Projected Change by 2040-2060 in Average Annual Precipitation per Wet Day (mm per day) at Eau Claire by Climate Projection 83  precipitation as increases. The NARCCAP projections of the annual amount of precipitation per wet day range from -0.5 to 1 mm across the different downscaling methods.  The projected change in annual frequency of wet days is largely dependent on downscaling method (Figure 3.7). Delta downscaling does not allow for changes in the frequency of wet days, as historical precipitation series are multiplied by a change factor (expressed as a ratio) to obtain future precipitation series. The NARCCAP simulations downscaled with the Raw approach project decreases in the frequency of wet days. The projections developed using the QM, MLR, and MLCR downscaling methods are approximately equally divided between increases and decreases in the frequency of wet days, with values ranging from -30 to 29 days per year. As described earlier, the delta climate projections were developed by applying change factors calculated by month to historical time series of temperature and precipitation. Thus, an understanding of the differences among the simulations in the magnitude of the monthly deltas is needed to interpret future changes in annual and seasonal means of the climate variables and the outputs of the yield model. This is particularly important for the temperature variables, given their contribution to freeze risk. The timing and magnitude of the monthly deltas of maximum and minimum temperature are shown for the CMIP5 models in Figures 3.8 and 3.9 respectively, and for the NARCCAP model combinations in Figure 3.10. There is large variability in the magnitude of the monthly temperature deltas and the timing for the largest deltas among the CMIP5 models. For example, compared to the other models, the MIROC-ESM and MIROC-ESM-CHEM models consistently project the largest increases in maximum and minimum temperature during late winter and early spring (i.e., February and March). Inter-model variability in the monthly magnitudes and fluctuations of the delta values is also seen for the 84  finer resolution NARCCAP model combinations, particularly for minimum temperatures, as seen by the large monthly deltas in January and February for the CRCM_CGCM3 and CRCM_CCSM simulations.   Figure 3.7: Projected Change by 2040-2060 in Average Annual Frequency of Wet Days at Eau Claire by Climate Projection 85    Figure 3.8: Monthly Deltas (°C) for Maximum Temperature from the CMIP5 Models The deltas were calculated for the model land gridpoint nearest the station location. The deltas were applied to observed series of maximum temperature to obtain future series of maximum temperature. The stations are shown by row, and the RCP by column. 86    Figure 3.9: Monthly Deltas (°C) for Minimum Temperature from the CMIP5 Models The deltas were calculated for the model land gridpoint nearest the station location. The deltas were applied to observed series of maximum temperature to obtain future series of maximum temperature. The stations are shown by row, and the RCP by column.  87     Figure 3.10: Temperature Deltas for the NARCCAP Simulations  Maximum and minimum temperature deltas (°C) per month used in downscaling the NARCCAP simulations using the delta approach. The stations are shown by row, with maximum temperature in the left column, and minimum temperature in the right column. 88  3.3.3 Projected Future Changes in Sour Cherry Production The projected future changes in sour cherry yield, buds remaining, poor pollination days, timing of critical phenological stages, and the frequency and severity of damaging freeze events, as simulated by the sour cherry model, are examined below for Eau Claire, Hart, and Maple City. Two measures of central tendency (the ensemble average and median) and two measures of dispersion (the full range of the projections and the range between the 25th and 75th percentiles) are used to summarize the projected changes and their uncertainty. 3.3.3.1 Changes in Sour Cherry Yield  Almost all of the projections developed from the CMIP5 models using the delta method, regardless of RCP, suggest a decrease in sour cherry yield at Eau Claire for 2040-2060 (Figure 3.11). Furthermore, all but one of the NARCCAP projections downscaled using the delta approach indicate a yield decrease. The uncertainty envelope is large, ranging for the CMIP5 projections from -41.9 to 7.1%, with an ensemble average (median) of -11.7 (-9.8)%, -9.0 (-6.3)%, and -13.4 (-12.9)% under RCPs 4.5, 6.0, and 8.5, respectively (Table 3.8). The yield decreases are much larger for the MIROC-ESM-CHEM and MIROC-ESM projections. As noted earlier, both of these models project large increases in maximum and minimum temperature in late winter and early spring. Removing the CMIP5 models with projected yield changes outside the 25th and 75th percentiles reduces the range greatly. The uncertainty range for the yield changes obtained from NARCCAP delta projections is less than half of that estimated from the CMIP5 projections, although the ensemble size is also smaller. However, the reduced 25-75th percentile range for the NARCCAP delta projections is within a few percent of that for the CMIP5 RCP 4.5 and 6.0 projections. WRFG_CCM3 is the only NARCCAP delta projection indicating a future increase in yield. This simulation had somewhat higher projections of  89     Figure 3.11: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Eau Claire by Climate Projection 90    Table 3.8: Simulated Changes (percent) by 2040-2060 in Sour Cherry Yield by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated yield between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles   Eau Claire CMIP5 Delta RCP 4.5 -11.7 -9.8 37.7 -10.8 17.0 CMIP5 Delta RCP 6.0 -9.0 -6.3 32.0 -7.1 9.4 CMIP5 Delta RCP 8.5 -13.4 -12.9 48.9 -12.2 18.4 NARCCAP Raw 2.3 1.3 35.8 0.6 11.0 NARCCAP Delta -7.2 -9.0 18.2 -8.1 11.0 NARCCAP QM 2.1 0.3 25.8 1.4 14.0 NARCCAP MLR 2.4 2.3 20.1 1.9 9.2 NARCCAP MLCR 3.6 0.0 31.1 1.4 17.5 Hart CMIP5 Delta RCP 4.5 -7.0 -5.8 37.0 -5.7 23.3 CMIP5 Delta RCP 6.0 -3.7 0.8 33.9 -1.4 15.9 CMIP5 Delta RCP 8.5 -10.5 -5.4 52.3 -7.6 26.8 NARCCAP Raw 2.5 -0.3 29.6 0.6 16.7 NARCCAP Delta 1.5 2.3 14.1 1.6 13.3 NARCCAP QM -1.9 -2.0 30.2 -1.6 26.5 NARCCAP MLR -0.6 -1.0 12.4 -1.0 4.9 NARCCAP MLCR 0.6 -3.1 33.6 -0.7 9.1 Maple City CMIP5 Delta RCP 4.5 -0.5 1.7 61.3 0.0 25.1 CMIP5 Delta RCP 6.0 -0.7 1.8 41.5 0.5 22.5 CMIP5 Delta RCP 8.5 -4.6 -2.1 58.8 -2.6 37.7 NARCCAP Raw 6.2 3.7 35.4 5.9 23.0 NARCCAP Delta 10.5 10.4 22.9 10.8 10.6 NARCCAP QM 14.7 14.5 36.3 14.4 34.1 NARCCAP MLR 0.9 0.1 15.4 0.0 4.2 NARCCAP MLCR 2.8 0.2 28.5 2.1 21.0     91  minimum temperature during spring and early summer for the land grid point nearest Eau Claire compared to the other NARCCAP simulations. Uncertainty in the sign of the projected change in yield at Eau Claire is much larger for the downscaling methods that, unlike the delta method, allow the variability of temperature to change in the future. The simulated future changes in yield when the NARCCAP simulations are downscaled by the Raw, QM, MLR, and MLCR methods are input into the cherry model are almost equally divided between increases and decreases in yield. Yield projections obtained from the more complex downscaling methods range from -13.5 to 25.8%, and the ensemble means are small although positive. Future yield projections obtained using the WRFG_CCSM simulation are consistently relatively large compared to the other NARCCAP simulations regardless of downscaling approach, but the type of downscaling applied appears to influence the relative magnitude of the estimated yield change for the other NARCCAP simulations.  The sign of the estimated yield changes at Hart obtained from the NARCCAP simulations downscaled using the raw, QM, MLR and MLCR methods is highly uncertain, similar to the findings for Eau Claire (Figure 3.12). Simulated yields range from -18.2 to 23.0%, and the ensemble means are close to zero (Table 3.8). However, unlike Eau Claire, the CMIP5 and NARCCAP projections downscaled using the delta method also suggest high uncertainty in the direction of the projected change. The future yields estimated from the CMIP5 models range from -42.7 to 9.6%, with an ensemble average (median) of -7.0 (-5.8)%, -3.7 (0.8)%, and -10.5 (-5.4)% for RCP 4.5, 6.0, and 8.5 respectively. The ensemble mean (median) of the yield changes obtained from the NARCCAP delta projections is 1.5 (2.3)%, although the uncertainty range is smaller compared to range obtained from the CMIP5 projections. As for Eau Claire, the projections obtained from the CMIP5 MIROC-ESM-CHEM and MIROC-ESM simulations are  92     Figure 3.12: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Hart by Climate Projection 93  associated with large yield decreases, and those downscaled from the NARCCAP WRFG_CCSM simulation are associated with relatively large yield increases. The yield estimates obtained from the NARCCAP simulations downscaled using the MLR and MLCR are particularly sensitive to outliers, as seen by the much smaller 25th - 75th percentile range compared to the full range.  As for Hart, approximately an equal number of the CMIP5 projections for Maple City lead to simulated yield decreases as increases (Figure 3.13), although the uncertainty range is larger (Table 3.8). The simulated yield changes obtained from CMIP5 projections range from -37.7 to 28.1%, and the ensemble mean (median) is -0.5 (1.7)%, -0.7 (1.8)%, and -4.6 (-2.1)% for RCP 4.5, 6.0, and 8.5, respectively. Even after removing the yield changes that fall outside the 25th and 75th percentiles, the uncertainty range for the CMIP5 projections is larger for Maple City compared to the other two locations. Unlike for Hart and Eau Claire, the pattern of the simulated yields obtained from the NARCCAP simulations downscaled using the delta methods differs from that of the yields obtained from the CMIP5 projections, with all but one of the NARCCAP delta projections associated with a yield increase. The majority of the simulated yield changes are also positive for the NARCCAP QM simulations, whereas the sign of the projected yield changes is inconsistent for the NARCCAP Raw, MLR and MLCR projections. Thus, the largest ensemble means are obtained from the NARCCAP Delta and NARCCAP QM projections.    94      Figure 3.13: Simulated Change (percent) by 2040-2060 in Average Sour Cherry Yield at Maple City by Climate Projection 95  3.3.3.2 Changes in Buds Remaining In the sour cherry model, yield is estimated from the buds remaining and the number of poor pollination days. Consequently, a strong association is expected between the estimated yield and estimated buds remaining, as seen for all 3 locations (Table 3.9). Interpretation of the sign of the future change in buds remaining is almost identical to that of estimated yield. At Eau Claire, almost all the climate projections downscaled using the delta method suggest a decrease in the average number of buds remaining by mid-century, whereas the sign of the change is mixed for the remaining climate projections (Figure 3.14). At Hart, the sign of the future change in buds remaining is highly uncertain, regardless of the type of climate projection used to estimate the change (Figure 3.15). The NARCCAP Delta and QM projections suggest an increase in buds remaining at Maple City by mid-century, whereas the number of positive and negative simulated changes is similar for the other downscaling methods (Figure 3.16). 3.3.3.3 Changes in Poor Pollination Days Differences between stations are once again highlighted by the simulated changes in poor pollination days. When the CMIP5 projections are input to the yield model, the majority of the projections, regardless of RCP, suggest an increase in the number of poor pollination days at Eau Claire (Figure 3.17). At Hart, an increase in the number in poor pollination days is simulated under RCPs 4.5 and 8.5, but not RCP 6.0 (Figure 3.18). This is in contrast to the mixed findings for simulated yield and buds remaining under all three RCPs, and also contrasts with increased number poor pollination days at Maple City for the RCP 6.0 climate projections and the uncertain sign of the projected changes for other two RCPs (Figure 3.19). The estimates of the number of poor pollination days from all the CMIP5 projections range from -0.3 to 2.5 days per year at Eau Claire, -0.5 to 3.3 days per year at Hart, and -0.8 to 2.2 days per year at Maple City 96  (Table 3.10). The sign of the projected change in the number of poor pollination days obtained from the NARCCAP simulations also varies by location. In particular, the majority of the estimated changes are positive for the NARCCAP delta projections at Eau Claire and Maple City but not Hart, whereas the NARCCAP QM simulations suggest a decrease in the number of poor pollination days at Eau Claire and Hart but an increase at Maple City.  Table 3.9: Simulated Changes (percent) by 2040-2060 in Buds Remaining by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated buds remaining between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 -12.8 -11.1 41.8 -12.4 18.5 CMIP5 Delta RCP 6.0 -10.7 -8.0 33.9 -9.3 9.3 CMIP5 Delta RCP 8.5 -14.9 -12.8 57.1 -13.2 20.2 NARCCAP Raw 1.9 0.2 35.1 -0.4 12.3 NARCCAP Delta -10.0 -10.6 21.8 -11.0 9.9 NARCCAP QM 2.1 -0.3 24.6 0.7 14.7 NARCCAP MLR 0.8 0.0 16.9 -0.3 3.2 NARCCAP MLCR 4.2 1.2 33.2 1.5 7.8 Hart CMIP5 Delta RCP 4.5 -7.5 -3.0 42.5 -4.9 22.7 CMIP5 Delta RCP 6.0 -3.7 -0.2 32.2 -1.2 12.8 CMIP5 Delta RCP 8.5 -12.0 -6.9 53.2 -8.8 29.3 NARCCAP Raw 2.5 0.1 29.7 1.7 24.9 NARCCAP Delta 1.4 2.1 15.5 1.6 11.1 NARCCAP QM -4.1 -1.8 25.5 -3.4 21.5 NARCCAP MLR 1.1 0.0 9.6 0.0 0.0 NARCCAP MLCR 1.6 0.6 29.1 0.1 10.6 Maple City CMIP5 Delta RCP 4.5 1.3 -0.3 70.4 2.9 32.6 CMIP5 Delta RCP 6.0 4.0 8.4 46.2 5.9 23.3 CMIP5 Delta RCP 8.5 -2.9 -4.4 72.6 -0.8 42.5 NARCCAP Raw 5.1 4.5 30.8 5.8 18.3 NARCCAP Delta 16.1 16.8 22.9 16.1 18.5 NARCCAP QM 18.0 21.4 54.4 18.7 27.5 NARCCAP MLR 2.3 0.4 18.9 0.3 2.0 NARCCAP MLCR 1.5 -2.8 26.1 0.1 22.4   97        Figure 3.14: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Eau Claire by Climate Projection  98    Figure 3.15: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Hart by Climate Projection 99         Figure 3.16: Simulated Change (percent) by 2040-2060 in Average Buds Remaining at Maple City by Climate Projection  100     Figure 3.17: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Eau Claire by Climate Projection 101       Figure 3.18: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Hart by Climate Projection  102         Figure 3.19: Simulated Changes by 2040-2060 in the Average Number of Poor Pollination Days per Year at Maple City by Climate Projection 103    Table 3.10: Simulated Changes by 2040-2060 in Poor Pollination Days by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated number of poor pollination days between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 1.1 1.1 2.5 1.0 1.0 CMIP5 Delta RCP 6.0 0.8 0.5 2.6 0.6 0.7 CMIP5 Delta RCP 8.5 1.2 1.3 2.2 1.2 0.9 NARCCAP Raw -0.2 -0.1 1.9 -0.2 1.5 NARCCAP Delta 0.5 0.5 1.2 0.5 1.0 NARCCAP QM -0.3 -0.5 2.0 -0.3 1.5 NARCCAP MLR -0.4 -0.5 2.6 -0.5 1.4 NARCCAP MLCR -0.1 0.1 2.9 0.0 2.1 Hart CMIP5 Delta RCP 4.5 0.8 0.8 2.9 0.8 1.2 CMIP5 Delta RCP 6.0 0.5 0.2 2.6 0.2 0.8 CMIP5 Delta RCP 8.5 1.0 1.1 3.6 0.9 1.8 NARCCAP Raw 0.0 0.2 2.3 0.1 1.1 NARCCAP Delta 0.3 0.2 1.6 0.3 1.0 NARCCAP QM -0.4 -0.5 2.1 -0.5 1.8 NARCCAP MLR 0.3 0.2 2.2 0.2 1.2 NARCCAP MLCR 0.1 0.1 1.5 0.1 1.4 Maple City CMIP5 Delta RCP 4.5 0.3 0.3 1.7 0.3 0.8 CMIP5 Delta RCP 6.0 0.5 0.6 1.8 0.6 0.5 CMIP5 Delta RCP 8.5 0.5 0.3 2.7 0.5 0.9 NARCCAP Raw -0.2 0.1 2.3 -0.2 2.0 NARCCAP Delta 0.2 0.2 0.5 0.2 0.3 NARCCAP QM -0.5 -0.5 2.0 -0.6 1.5 NARCCAP MLR 0.1 0.1 2.6 0.0 1.8 NARCCAP MLCR -0.4 -0.7 2.6 -0.4 2.0     104  3.3.3.4 Changes in Phenology The timing of phenology plays an important role in freeze risk and thus influences the simulated yields. All the climate projections point to an earlier occurrence of phenological stage 2 (side green) (Table 3.11) at Eau Claire (Figure 3.20), Hart (Figure 3.21), and Maple City (Figure 3.22), regardless of the simulation source, downscaling method, or greenhouse gas emissions. For a particular projection type, the ensemble means of the simulated shift in the timing of stage 2 are nearly identical for the three locations, although differences in the ensemble means are seen between projection types. The simulated changes are somewhat larger for the CMIP5 RCP 4.5 and 8.5 projections, with the ensemble means varying across the stations from -16.8 to -15.2 days for RCP 4.5 and from -20.9 to -18.6 days for RCP 6.0. The ensemble means for the other projection types are smaller, ranging from -12.9 to -9.0 days. At all stations, the simulated advancement of stage 2 is much larger for the climate projections developed from the CMIP5 MIROC-ESM-CHEM and MIROC-ESM models, highlighting the influence of the choice of GCM models to include in an ensemble on the ensemble range. For the NARCCAP-derived climate projections, those derived from the WRFG-CCSM, HRM3_GFDL, and CRCM_CGCM3 are more likely to be associated with a larger simulated advancement in stage 2.  At all stations, the simulated future changes in the time of occurrence of phenological stage 8 (full bloom) are similar to those for stage 2, suggesting a uniform advancement of phenological stages by the mid-century across the study region (Table 3.12 and Figures 3.23-3.25). Again, the range of the projected changes at an individual station is large, but the differences between stations are small. Removing the projections that fall outside the 25th-75th percentiles substantially reduces the range, especially for the CMIP5 projections, but has only a 105  small influence on the ensemble means. The ensemble means for the advancement of stage 8 vary across the three stations from -15.2 to -14.0 days for CMIP5 RCP 4.5, from -18.6 days to -17.3 days for CMIP5 8.5, and from  13.2 to -9.1 days for the remaining climate projections. The smallest estimated changes (-10.1 to -9.7 days) are obtained from the NARCCAP MLR and MLCR projections. Table 3.11: Simulated Changes by 2040-2060 in Date of Phenological Stage 2 (Side Green) by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated dates between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 -15.4 -14.8 25.2 -14.4 11.5 CMIP5 Delta RCP 6.0 -11.8 -10.3 23.0 -9.8 6.8 CMIP5 Delta RCP 8.5 -18.6 -16.0 35.0 -16.6 13.1 NARCCAP Raw -9.4 -10.1 9.5 -9.9 6.1 NARCCAP Delta -11.5 -12.2 9.9 -11.8 6.2 NARCCAP QM -10.8 -11.0 8.2 -11.1 5.2 NARCCAP MLR -8.2 -8.0 7.6 -8.2 6.9 NARCCAP MLCR -8.9 -9.4 9.4 -9.1 2.7 Hart CMIP5 Delta RCP 4.5 -15.2 -14.2 23.5 -14.2 10.8 CMIP5 Delta RCP 6.0 -11.7 -8.7 22.9 -9.3 6.9 CMIP5 Delta RCP 8.5 -19.2 -17.7 37.5 -17.3 14.7 NARCCAP Raw -10.1 -11.1 9.4 -10.8 3.3 NARCCAP Delta -10.3 -10.5 11.3 -10.6 6.9 NARCCAP QM -11.6 -12.5 7.6 -11.9 5.0 NARCCAP MLR -9.4 -9.0 7.3 -9.6 4.2 NARCCAP MLCR -9.8 -9.7 7.0 -9.9 5.3 Maple City CMIP5 Delta RCP 4.5 -16.8 -15.6 35.3 -15.0 10.2 CMIP5 Delta RCP 6.0 -12.9 -9.0 29.9 -10.2 9.7 CMIP5 Delta RCP 8.5 -20.9 -18.5 45.5 -18.3 15.5 NARCCAP Raw -9.3 -10.3 10.1 -9.8 4.8 NARCCAP Delta -10.4 -11.4 10.9 -10.9 6.0 NARCCAP QM -10.4 -9.9 8.7 -10.4 4.7 NARCCAP MLR -8.9 -8.7 6.9 -9.0 6.1 NARCCAP MLCR -9.7 -9.5 7.2 -9.7 5.9   106       Figure 3.20: Simulated Change by 2040-2060 in Average Date of Stage 2 at Eau Claire by Climate Projection 107      Figure 3.21: Simulated Change by 2040-2060 in Average Date of Stage 2 at Hart by Climate Projection 108      Figure 3.22: Simulated Change by 2040-2060 in Average Date of Stage 2 at Maple City by Climate Projection 109    Table 3.12: Simulated Changes by 2040-2060 in Date of Phenological Stage 8 (Full Bloom) by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated dates between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 -14.6 -12.7 22.5 -13.5 10.3 CMIP5 Delta RCP 6.0 -11.3 -9.3 22.0 -9.1 6.7 CMIP5 Delta RCP 8.5 -17.7 -16.5 34.1 -15.9 9.6 NARCCAP Raw -9.8 -10.0 8.7 -10.1 4.8 NARCCAP Delta -10.2 -11.0 10.3 -10.7 5.2 NARCCAP QM -11.5 -11.9 9.5 -11.8 6.4 NARCCAP MLR -9.1 -9.2 9.3 -9.3 3.3 NARCCAP MLCR -9.6 -9.7 6.3 -9.7 4.5 Hart CMIP5 Delta RCP 4.5 -14.0 -12.8 19.0 -12.9 7.0 CMIP5 Delta RCP 6.0 -10.9 -8.9 18.5 -8.9 5.0 CMIP5 Delta RCP 8.5 -17.3 -15.5 31.6 -15.3 9.7 NARCCAP Raw -10.7 -11.1 9.4 -10.8 6.3 NARCCAP Delta -11.0 -11.6 10.4 -11.5 5.7 NARCCAP QM -13.2 -15.3 13.7 -13.8 9.5 NARCCAP MLR -9.9 -9.9 8.3 -10.0 5.0 NARCCAP MLCR -10.1 -10.2 9.4 -10.6 3.7 Maple City CMIP5 Delta RCP 4.5 -15.2 -14.0 31.7 -13.1 6.4 CMIP5 Delta RCP 6.0 -11.8 -7.8 27.8 -8.9 6.8 CMIP5 Delta RCP 8.5 -18.6 -15.9 41.6 -15.5 11.9 NARCCAP Raw -10.2 -11.0 9.4 -10.5 7.0 NARCCAP Delta -10.5 -11.5 8.1 -11.1 4.5 NARCCAP QM -11.5 -13.0 11.3 -12.0 6.4 NARCCAP MLR -9.7 -9.7 6.7 -9.9 4.0 NARCCAP MLCR -9.9 -10.5 8.9 -10.4 3.9      110        Figure 3.23: Simulated Change by 2040-2060 in Average Date of Stage 8 at Eau Claire by Climate Projection  111      Figure 3.24: Simulated Change by 2040-2060 in Average Date of Stage 8 at Hart by Climate Projection 112      Figure 3.25: Simulated Change by 2040-2060 in Average Date of Stage 8 at Maple City by Climate Projection 113  3.3.3.5 Changes in Frequency and Severity of Damaging Freeze Events The frequency and severity of damaging freeze events, as simulated by the sour cherry model, are also examined. A damaging freeze is defined as any freeze event that results in bud loss. The severity of a damaging freeze is the average damage amount per damaging freeze event, with a value of 0 indicating no loss of buds and a value of 1 indicating complete bud loss. The uniformity of the sign of the projected changes in phenology and the small differences between stations highlight that the uncertainty and inter-station differences discussed earlier for simulated buds remaining, and subsequently sour cherry yield, are primarily a function of differences between climate projections and locations in the simulated frequency and severity of damaging freezes. The ensemble plots of the frequency of damaging freezes for Eau Claire, Hart, and Maple City support this interpretation. At Eau Claire, the simulated changes in the frequency of damaging freezes are positive for all but a few (2 or 3 depending on RCP) delta-downscaled CMIP5 projections, with the increases as large as 80 days per year for some of the RCP 8.5 projections (Figure 3.26). The simulated increases in the frequency of damaging freezes tend to be larger for those CMIP5 projections with the largest advancement in phenological stage (e.g., MIRO-ESM, MIROC-ESM-CHEM). Also, all but one of the delta-downscaled NARCCAP simulations project an increase in the frequency of damaging freezes. On the other hand, the simulated frequencies of damaging freeze events obtained from climate projections derived from the NARCCAP simulations using downscaling methods that allow the variability of temperature to change in the future suggest either little change in the frequency of damaging freezes or that a decrease in frequency is as likely as an increase. These differences are clearly seen in the ensemble means which are positive for the delta-downscaled projections, but are close to zero for the other projections (Table 3.13). 114     Figure 3.26: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Eau Claire by Climate Projection 115    Table 3.13: Simulated Changes by 2040-2060 in Frequency of Damaging Freezes by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated frequency between the 25th and 75th percentiles are also shown. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 15.7 9.5 42.0 15.0 30.0 CMIP5 Delta RCP 6.0 10.1 9.0 41 8.9 16 CMIP5 Delta RCP 8.5 21.9 12.0 104 16.3 31 NARCCAP Raw 1.0 -1.0 23 0.5 16 NARCCAP Delta 11.8 12.0 29 12.2 13 NARCCAP QM 0.8 2.5 28 1.3 21 NARCCAP MLR -4.5 0.0 47 0.2 3 NARCCAP MLCR -7.5 0.5 72 -0.7 6 Hart CMIP5 Delta RCP 4.5 10.9 5.5 52 7.9 26 CMIP5 Delta RCP 6.0 4.8 1.5 32 2.3 13 CMIP5 Delta RCP 8.5 17.3 6.5 71 11.4 32 NARCCAP Raw 2.4 0.5 9 2.0 9 NARCCAP Delta -1.5 -3.0 17 -1.8 14 NARCCAP QM 4.8 3.0 32 3.3 22 NARCCAP MLR -3.0 0.0 26 0.0 0 NARCCAP MLCR -1.4 0.0 22 -0.5 9 Maple City CMIP5 Delta RCP 4.5 2.1 1.0 87 -0.8 32 CMIP5 Delta RCP 6.0 -4.1 -8.5 65 -6.4 21 CMIP5 Delta RCP 8.5 5.3 6.5 86 1.2 47 NARCCAP Raw 1.4 0.5 20 0.8 2 NARCCAP Delta -18.4 -20.5 19 -18.7 16 NARCCAP QM -16.9 -16.5 37 -17.0 25 NARCCAP MLR -2.8 0.0 22 -0.3 3 NARCCAP MLCR -1.3 2.0 29 -0.5 21     116  In contrast, all the different types of climate projections suggest that the direction of the change in the frequency of damaging freezes is uncertain at Hart (Figure 3.27). The ensemble means for the CMIP5 projections are positive, even though an approximately equal number of models project an increase as a decrease in damaging freeze frequency, as the magnitudes of the projected changes are larger for those models with an increase in frequency (Table 3.13). The ensemble means for the projections derived from the NARCCAP simulations are close to zero. Maple City is the only location for which some of the climate projections suggest a substantial future decrease in the frequency of damaging freezes (Figure 3.28). For example, almost all of the projections derived from the NARCCAP simulations using the Delta and QM downscaling methods suggest that damaging freezes will decrease in the future. For the other NARCCAP-derived projections, the sign of the projected change is uncertain. In terms of the CMIP5 delta projections, a majority of the projections derived from the CMIP5 RCP 6.0 simulations suggest a decrease in the frequency of damaging freezes, whereas the direction of the change is unclear for the CMIP5 RCP 4.5 and 8.5 projections. However, in contrast to Hart, projected increases and decreases are more similar in magnitude, resulting in smaller ensemble means for the CMIP5 projections for Maple City (Table 3.13).  The ensemble plots of changes in the average severity of individual damaging freezes for the CMIP5 projections suggest little change in severity at Eau Claire (Figure 3.29), somewhat larger changes at Hart although the direction of the change is unclear (Figure 3.30), and the potential for substantially greater severity at Maple City (Figure 3.31). With the exception of the NARCCAP delta projections, the interpretation of future changes in the severity damaging freeze events is not as straightforward for the NARCCAP-derived projections, with considerable variability across the projections in the magnitude and direction of the simulated change for all 117  three stations (Table 3.14). Some NARCCAP simulations downscaled using the Raw and MLR methods do not project any damaging freezes in the mid-century future time slice.     Figure 3.27: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Hart by Climate Projection 118      Figure 3.28: Simulated Change by 2040-2060 in Frequency (Number of Days per Year) of Damaging Freezes at Maple City by Climate Projection 119      Figure 3.29: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Eau Claire by Climate Projection The climate models with no projected change shown in the figure did not project any damaging freeze events during the 21-year period, 2040-2060.  120        Figure 3.30: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Hart by Climate Projection The climate models with no projected change shown in the figure did not project any damaging freeze events during the 21-year period, 2040-2060.  121      Figure 3.31: Simulated Change (percent) by 2040-2060 in Average Severity of Damaging Freezes at Maple City by Climate Projection The climate models with no projected change shown in the figure did not project any damaging freeze events during the 21-year period, 2040-2060. 122    Table 3.14: Simulated Changes (percent) by 2040-2060 in Average Severity of Damaging Freezes by Climate Projection Type Ensemble average, median, and range by source of climate simulation and downscaling method are provided. The ensemble average and range of the estimated severity between the 25th and 75th percentiles are also shown. Some climate models projected few or no damaging freeze events during the 21-year period, affecting the calculations of average severity of damaging freezes. Station Downscaling Method Ensemble Average Ensemble Median Ensemble Range Ensemble Average for 25-75th Percentiles Ensemble Range for 25-75th Percentiles Eau Claire CMIP5 Delta RCP 4.5 5.7 6.8 25.7 5.8 7.8 CMIP5 Delta RCP 6.0 6.3 7.1 21.2 6.8 6.7 CMIP5 Delta RCP 8.5 2.8 4.4 39.1 3.1 12.5 NARCCAP Raw 26.9 22.9 118.3 26.9 43.4 NARCCAP Delta 6.9 6.7 19.4 6.9 12.2 NARCCAP QM -0.9 -4.8 83.4 -6.1 34.2 NARCCAP MLR -62.2 -100.0 151.2 -62.2 0.0 NARCCAP MLCR -3.9 -5.5 48.5 -4.8 32.7 Hart CMIP5 Delta RCP 4.5 5.1 4.4 44.5 5.4 23.0 CMIP5 Delta RCP 6.0 -0.5 0.6 37.7 -1.2 20.0 CMIP5 Delta RCP 8.5 10.7 9.1 60.3 10.8 25.0 NARCCAP Raw -47.1 -55.1 121.8 -47.1 89.7 NARCCAP Delta 15.1 16.8 25.5 15.5 14.5 NARCCAP QM 4.5 2.4 43.3 2.3 16.2 NARCCAP MLR -100.0 -100.0 - - - NARCCAP MLCR 4.2 3.3 69.8 3.3 38.0 Maple City CMIP5 Delta RCP 4.5 13.4 17.7 49.4 14.8 27.7 CMIP5 Delta RCP 6.0 6.1 8.0 32.3 5.6 17.9 CMIP5 Delta RCP 8.5 18.1 17.1 60.4 19.1 17.3 NARCCAP Raw -2.2 3.6 68.3 -2.2 7.2 NARCCAP Delta 7.4 8.6 14.1 8.2 7.3 NARCCAP QM -2.0 -5.7 32.8 -1.8 21.4 NARCCAP MLR -8.1 -30.5 215.6 -8.1 52.8 NARCCAP MLCR 13.4 9.6 98.9 8.7 31.7    123  3.4 Summary and Discussion This study examined the historical trends and projected future changes in sour cherry yield and several related variables (i.e., timing of sensitive phenology stages, frequency and severity of damaging freezes, the number of poor pollination days, and the number of buds remaining), as simulated by a physically-based empirical yield model for three locations in western Michigan (Eau Claire, Hart, and Maple City). The study employed a large ensemble of GCM simulations from the CMIP5 archive and RCM simulations from the NARCCAP archive. Both sets of simulations were downscaled to the location of the three stations. The major findings of the study are as follows. 1)  Linear trends in simulated sour cherry yield and related variables are mostly insignificant for the historical period of 1960-2015. Exceptions occur at Eau Claire in southwestern Michigan, where a significant positive trend in the frequency of damaging freeze events and a significant advancement in the timing of phenological stage 2 (side green) and stage 8 (full bloom) were found. Schultze et al. (2014) earlier found little change or even a slight decrease in the frequency of late frosts during the grape growing season in southwestern Michigan. Therefore, the increases in simulated frequency of damaging freezes at Eau Clare are likely the result of the significant advancement in the timing of sensitive phenological stages and the occurrence of freezes during these sensitive stages. This is in agreement with the finding by Augspurger (2013) that frost damage to various woody species in neighboring Illinois increased since 1980 due to advancing phenology. There is some evidence of nonlinear trends and/or breakpoints in the simulated time series. At all three stations, the simulated number of poor pollination days and damaging freezes increased during the last decade. Also, inspection of the time series of simulated yield for a longer period of record, extending back to 1894, at Hart 124  suggests periods of the record were dominated by increases in yield whereas other periods were dominated by negative yield trends.  2) The large of ensemble of climate projections suggests that annual maximum and minimum temperature in western Michigan will increase by the mid-21st century, with projected changes ranging from 1 to 4.6°C for the projections derived from the CMIP5 simulations and from 1.2 to 3.2°C for the projections derived from the NARCCAP simulations. The majority of the downscaled CMIP5 and NARCCAP simulations project an increase in the amount of precipitation per wet day. The sign of the projected change in the frequency of wet days is inconsistent among the downscaled NARCCAP simulations.  3) The timing of critical phenological stages is expected to advance to earlier dates by the mid-21st century. Although the amount of change at any particular location varies with the type of climate projection, the ensemble means are similar for the three locations and two phenological stages, suggesting a uniform advancement of phenology across the study area. This result is consistent with other studies which projected earlier dates of phenology in Europe (Chmielewski et al., 2010; Eccel et al., 2009; Hoffmann and Rath, 2013; Molitor et al., 2014).  4) The sign of the projected changes in the frequency and severity of damaging freezes is inconsistent, varying by climate projection type and location. In general, the climate projections downscaled from the NARCCAP simulations using methods that allow temperature variability to change in the future suggest little change in the frequency of damaging freezes or that a decrease in frequency is as likely as an increase, regardless of location. The CMIP5 and NARCCAP simulations downscaled using the simple delta method that assumes constant variance vary by location. The frequency of damaging freezes is expected to increase at Eau Claire, whereas the direction of the projected change is uncertain at Hart and Maple City. However, the average 125  severity of individual damaging freezes is more likely to increase at Maple City compared to the other two locations. These findings contrast with those of Rochette et al. (2004) who found that advancement in the phenology of fruit trees in southern Ontario led to projected increases in frost damage risk. Rather, the results of this study indicate that the relationship between the timing of sensitive phenological stages and the frequency of cold temperatures is complex and varies spatially. This interpretation supports that of Winkler et al. (2012b) who concluded, based on a large ensemble of downscaled CMIP3 simulations, that future changes in frost risk for sour cherry production in Michigan are uncertain.  5) A unique contribution of the study is consideration of future changes in the number of poor pollination days. As for the frequency of damaging freezes, the ensemble of climate projections indicates considerable uncertainty in future changes, although, in general, the frequency of poor pollination days appears more likely to increase in the future at Eau Claire and Maple City compared to Hart.  6) The ensemble of climate projections used in this study highlights the uncertainty surrounding projections in future sour cherry yield. This finding is in concordance with the large uncertainty in the projected frequency of damaging freeze events. The focus on the frequency and severity of damaging freeze events on yield is unique to this study. Previous analyses, such as that of Rochette et al. (2004), for the most part estimated freeze risk solely on the timing of the last freeze relative to phenological development. This approach is overly simplistic as the buds remaining on sour cherry trees at the time of last spring frost, and the consequent impact on yield, represents the cumulative impact of multiple freeze events during the phenological development of the buds.   126  7) The use of a large and varied ensemble of climate projections is essential to better understand the uncertainty associated with projected future changes. Also, ensemble means can be misleading and need to be accompanied by depictions of the range of the projected changes by the ensemble members. Furthermore, the ensemble needs to include projections developed from different GCM or RCM simulations, greenhouse gas concentrations, and downscaling methods. Considerable differences were observed between the projections downscaled using the simple delta approach that do not allow for changes in variability of precipitation frequency and the dynamically-downscaled projections that were further downscaled using more complex downscaling approaches. Also, the uncertainty envelope was generally larger for the projections downscaled from the CMIP5 GCM simulations than for those downscaled from the NARCCAP RCM simulations. This finding agrees with that of Kunkel et al. (2013) and Pryor et al. (2013), both of whom found a larger range in projected future temperature change from CMIP3 GCM simulations than from the NARCCAP RCM simulations. The development and use of a large ensemble of climate projections contrasts with a number of earlier studies. Some studies of future changes in perennial crop production used a single climate projection (e.g., Chmielewski et al. 2010; Kaukoranta et al., 2010; Ladányi et al., 2010), and those that did include multiple projections often only accounted for one source of uncertainty in the climate projections such as different greenhouse gas emissions (e.g., Eccel et al., 2009; Rochette et al., 2004). Few studies considered uncertainty introduced by the choice of downscaling method. This study highlights that climate projections can vary substantially depending on emissions scenario, downscaling method, and climate model, indicating that more inclusive climate projection ensembles should be used more often in climate impact studies.  127  As with any study, there are a number of limitations, but one aspect that particularly requires further attention is the parameter and structural uncertainty of the yield model. The model uses temperature thresholds derived from a chamber study to estimate phenological stage and sensitivity to freeze damage, but chamber conditions are not necessarily representative of conditions in orchards. Furthermore, the relationships estimated between temperature and phenology, as well as between pollination, freezes, and yield, from historical conditions are assumed to hold true in the future. In other words, the model relationships are assumed to be stationary, which may not be the case. Spatial variations in model performance were also observed. The simulated yield at Maple City explained 31% of the variation of observed regional yields in northwest Michigan. However, the yield model performed less well at Hart as compared to observed regional yields for west central Michigan. The poor performance at this station may be due to microclimate influences at Hart, causing the station to be unrepresentative of the region. Additionally, observed yields are unavailable for southwestern Michigan to validate the model at Eau Claire. Another complication is that, because of microclimate influences, the three stations selected to represent the fruit production areas in western Michigan based on the length and quality of their climate records, may not be representative of the majority of orchard locations.  While there are many sources of uncertainty in this study, the implications of the projected future changes are important for fruit growers in Michigan. Growers need to be aware of the uncertainty surrounding future projections so that the decisions they make are robust and flexible.   128  CHAPTER 4: CONCLUSION This research focused on the influence of climate on sour cherry trees in the Great Lakes Region, specifically the negative impacts of springtime freezes on production. The historical relationship was examined in Chapter 2 and the historical and future temporal trends were analyzed in Chapter 3.  In Chapter 2 the timing and physical characteristics of springtime freezes were examined for 1960-2015 with respect to sour cherry production in the Great Lakes Region. A sour cherry model was used to simulate phenological growth, and temperature thresholds were used to estimate freeze damage to tree buds. It was found that in general springtime freezes were most common in early spring and in the north. Classic radiation freezes were more common in later spring, and overall only a slight majority of the freezes that were estimated to cause damage at the stations in the northwest and west central regions were found to be classic radiation events with calm, clear conditions at night. At the southwest station the majority of the simulated damaging freezes were found to be non-radiation freezes, possibly due microclimate effects. Most of the damaging freeze events were found to occur towards the end of the spring season because fruit trees are more vulnerable to freeze damage at later stages of phenology. Data from reference stations showed the northwest region experienced the highest frequency of damaging freezes compared to the west central and southwest regions. Moreover, data from additional stations highlighted a spatial trend of more severe springtime freezes in the north and inland from Lake Michigan. These results highlight the importance of microclimate effects and consequently the significance of orchard location. Furthermore, the characteristics of simulated damaging freezes were examined, which is unique to this study. The results showed that 129  damaging freezes were more common in the later phenological stages, but the freezes during the earlier stages caused more damage on average.   The sour cherry model was evaluated by comparing the results of station climatic data series with regional yields for northwest and west central Michigan. The sour cherry model at the station Maple City accounted for 30 percent of the variation of the observed yields in the northwest region. The results at the station Hart did not perform well against the west central observed yields. Because there is a lack of long series of historical climatic data, gridded datasets were used in the sour cherry model to determine if they would be suitable for this type of application. It was found that NLDAS gridded data did not account for any of the interannual variability in yields and was determined to be unsuitable for this application. The gridded dataset PRISM accounted for some of the variability in yields, but additional point data would be the most helpful in this type of application. In Chapter 3, a yield algorithm based on springtime freeze damage and pollination conditions was used to examine temporal trends in sour cherry yield. Historical trends in simulated yield and the underlying factors, springtime freeze damage, pollination environment, and phenology, were analyzed. Furthermore, an ensemble of climate projections was used to determine how yield may change in the future, by the mid-century time slice 2040-2060. This chapter focused on the historical trends and projected changes in yield at three stations, Maple City, Hart, and Eau Claire, which span from north to south respectively, along the west coast in Michigan.    The linear historical temporal trends were generally insignificant for the simulated variables, except for phenology at Eau Claire which indicated a statistically significant 130  advancement of dates in phenological stage 2 (side green) and stage 8 (full bloom). The temporal trends in simulated yield and related variables showed some non-linearity in the series.   To examine changes under the projected future climate, a large ensemble was used with 88 individual members to account for uncertainty in the projections. The climate scenarios projected an increase in maximum and minimum temperature for all stations. Due to the projected increases in temperature, phenology was also projected to advance. However, there was a wide range of future outcomes in changes in simulated sour cherry yield, due to varying changes in frequency of damaging freezes and poor pollination days by climate projection.   The large ensemble of climate scenarios highlighted the sources of uncertainty in climate projections, through RCP or emissions scenario, downscaling method, and climate model. All these factors played a large role in the resulting projections of changes in sour cherry yield at each of the stations. The scenarios under RCP 8.5, which represents the highest radiative forcing by the end of the century, projected the most extreme changes in the variables, compared to RCP 4.5 and 6.0. The different downscaling methods allowed for varying amounts changes in climate variability in the future. The projected changes varied by scenarios downscaled with Delta, Raw, QM, MLR, and MLCR methods. Generally, the CMIP5 GCMs projected a larger of changes than the NARCCAP RCMs, and removing the more extreme scenarios did influence the ensemble averages and ranges.   There are some limitations to this research, which are discussed in more detail in the individual chapters. Due to lack of long climatic series, only three stations were used to represent the Great Lakes production areas. Furthermore, these stations may not accurately represent the weather at orchards, as microclimate has a strong influence in the region. This limitation was highlighted in the validation of the sour cherry model, as it did not perform well at Hart.  131   The results from both Chapter 2 and 3 have implications for fruit growers in the Great Lakes Region, as well as valuable information regarding climate impact assessments. Results from Chapter 2 showed that there are not as es are previously thought. It is important for growers to be aware that damaging freezes are most common during the later phenological stages, due to the increasing vulnerability of the fruit crops. However, when freezes occur in the early phenological stages, the damage is generally more severe due to colder temperatures. Overall, it was found that microclimate played a strong role in the temperatures during freezes, so orchard location is very important in preventing freeze damage. Furthermore, growers should consider robust adaptation measures to climate change as this study showed that there is a large amount of uncertainty surrounding how yield will change in the future. This results of this study indicated that the use of gridded data in climate assessments should be carefully considered as it may not be suitable in regions where microclimate is important. Moreover, this study highlighted that there is a large range in uncertainty in future climate projections, indicating the importance of using a large, inclusive ensemble of climate projections when studying the impacts of the future climate.    132             APPENDIX   133  Cropscape  The CropScape data layers for Michigan in 2013 were developed from satellite imagery from the Landsat 8 OLI/TIRS sensor, Landsat 7 ETM+ sensor, and the Disaster Monitoring Constellation (DMC) DEIMOS-1 and UK2 sensors. Additional information from the United States Geological Survey (USGS) National Elevation Dataset (NED) and the imperviousness and canopy data layers from the USGS National Land Cover Database (NLCD) in 2006 were used in the CropScape data layers. The data layers were validated by the USDA Farm Service Agency (FSA) Common Land Unit (CLU) Program data. The metadata states that the 2013 Michigan classification matches the ground truth data, the omissions error is when a pixel is excluded from the category to which it belongs, and the commission error is when a pixel is included in an incorrect category.    134             BIBLIOGRAPHY   135  BIBLIOGRAPHY  Abraham, Z., Tan, P.-N., Perdinan, Winkler, J.A., Zhong, S., Liszewska, M., 2014. Contour  regression: A distribution-regularized regression framework for climate modeling. Statistical Analysis and Data Mining 7, 272281. doi:10.1002/sam.11222  Alexandersson, H., 1986. A homogeneity test applied to precipitation data. Journal of  climatology 6, 661675.  Andresen, J.A., Alagarswamy, G., Ritchie, J.T., Rotz, C.A., LeBaron, A.W., 2001. Assessment  of the impact of weather on maize, soybean, and alfalfa production in the Upper Great Lakes Region of the United States, 18951996. Agronomy Journal 93, 10591070. doi:10.2134/agronj2001.9351059x  Andresen, J., Hilberg, S., Kunkel, K., 2012. Andresen, J., S. Hilberg, K. Kunkel, 2012: Historical  Climate and Climate Trends in the Midwestern USA. In: U.S. National Climate Assessment Midwest Technical Input Report. J. Winkler, J. Andresen, J. Hatfield, D. Bidwell, and D. Brown, coordinators. Available from the Great Lakes Integrated Sciences and Assessments (GLISA) Center, http://glisa.msu.edu/docs/NCA/MTIT_Historical.pdf., U.S. National Climate Assessment Midwest Technical Input Report. Great Lakes Integrated Sciences and Assessments (GLISA) Center.  Augspurger, C.K., 2013. Reconstructing patterns of temperature, phenology, and frost damage  over 124 years: Spring damage risk is increasing. Ecology 94, 4150. doi:10.1890/12-0200.1  Barden, J.A., Neilsen, G.H., 2003. Selecting the orchard site, site preparation and orchard  planning and establishment, in: Apples: Botany, Production and Uses. CABI Publishing, Wallingford, UK, pp. 237265.  Buishand, T.A., 1982. Some methods for testing the homogeneity of rainfall records. Journal of  hydrology 58, 1127.  CCSP, 2008: Climate Models: An Assessment of Strengths and Limitations. A Report by the  U.S. Climate Change Science Program and the Subcommittee on Global Change Research [Bader D.C., C. Covey, W.J. Gutowski Jr., I.M. Held, K.E. Kunkel, R.L. Miller, R.T. Tokmakian and M.H. Zhang (Authors)]. Department of Energy, Office of Biological and Environmental Research, Washington, D.C., USA, 124 pp.  Charrier, G., Ngao, J., Saudreau, M., Améglio, T., 2015. Effects of environmental factors and  management practices on microclimate, winter physiology, and frost resistance in trees. Front. Plant Sci. 259. doi:10.3389/fpls.2015.00259  136  Chmielewski, F.-M., 2013. Phenology in Agriculture and Horticulture, in: Schwartz, M.D. (Ed.),  Phenology: An Integrative Environmental Science. Springer Netherlands, pp. 539561.  Chmielewski, F.-M., Blümel, K., Henniges, Y., 2010. Climate Change and Possible Late Frost  Damages to Apple Trees in Germany, in: Berichte Des Meteorologischen Instituts Der Albert-Ludwigs-Universität Freiburg. Presented at the Proceedings of the 7th Conference on Biometeorology, p. 50.  Chmielewski, F.-M., Müller, A., Bruns, E., 2004. Climate changes and trends in phenology of  fruit trees and field crops in Germany, 19612000. Agricultural and Forest Meteorology 121, 6978. doi:10.1016/S0168-1923(03)00161-8  Chuine, I., de Cortazar-Atauri, I.G., Kramer, K., Hänninen, H., 2013. Plant Development  Models, in: Schwartz, M.D. (Ed.), Phenology: An Integrative Environmental Science. Springer Netherlands, Dordrecht, pp. 275293.  Cittadini, E.D., de Ridder, N., Peri, P.L., van Keulen, H., 2006. A method for assessing frost  damage risk in sweet cherry orchards of South Patagonia. Agricultural and forest meteorology 141, 235243. doi:10.1016/j.agrformet.2006.10.011  Dai, J., Wang, H., Ge, Q., 2013. The decreasing spring frost risks during the flowering period for  woody plants in temperate area of eastern China over past 50 years. J. Geogr. Sci. 23, 641652. doi:10.1007/s11442-013-1034-6  Daly, C., Halbleib, M., Smith, J.I., Gibson, W.P., Doggett, M.K., Taylor, G.H., Curtis, J.,  Pasteris, P.P., 2008. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. International journal of climatology 28, 20312064. doi:10.1002/joc.1688  Daly, C., Widrlechner, M.P., Halbleib, M.D., Smith, J.I., Gibson, W.P., 2012. Development of a  New USDA Plant Hardiness Zone Map for the United States. Journal of Applied Meteorology and Climatology 51, 242. doi:10.1175/2010JAMC2536.1  Dennis, F.G., Howell, G.S., 1974. Cold Hardiness of Sour Cherry Bark and Flower Buds (No.  220). Michigan State University, Agricultural Experiment Station.  Eccel, E., Rea, R., Caffarra, A., Crisci, A., 2009. Risk of spring frost to apple production under  future climate scenarios: the role of phenological acclimation. International journal of biometeorology 53, 273286. doi:10.1007/s00484-009-0213-8  Fitchett, J.M., Grab, S.W., Thompson, D.I., 2015. Plant phenology and climate change Progress  in methodological approaches and application. Progress in Physical Geography 39, 460482. doi:10.1177/0309133315578940  Fitchett, J.M., Grab, S.W., Thompson, D.I., Roshan, G., 2014. Increasing frost risk associated  137  with advanced citrus flowering dates in Kerman and Shiraz, Iran: 19602010. International journal of biometeorology 58, 18111815. doi:10.1007/s00484-013-0778-0  Gombos, B., Koles, P., Montvajszki, M., 2011. Spatial differences of night temperature in hilly  regions and its horticultural importance. Agriculture and Environment 3, 102109.  Goyette, S., McFarlane, N.A., Flato, G.M., 2000. Application of the Canadian regional climate  model to the Laurentian great lakes region: Implementation of a lake model. Atmosphere-Ocean 38, 481503. doi:10.1080/07055900.2000.9649657  Gu, L., Hanson, P.J., Post, W.M., Kaiser, D.P., Yang, B., Nemani, R., Pallardy, S.G., Meyers, T.,  2008. The 2007 Eastern US Spring Freeze: Increased Cold Damage in a Warming World. BioScience 58, 253262. doi:10.1641/B580311  Guédon, Y., Legave, J.M., 2008. Analyzing the time-course variation of apple and pear tree dates  of flowering stages in the global warming context. Ecological Modelling 219, 189199. doi:10.1016/j.ecolmodel.2008.08.010  Hanssen-Bauer, I., Achberger, C., Benestad, R., Chen, D., Førland, E., Faculty of Sciences,  Naturvetenskapliga fakulteten, Department of Earth Sciences, University of Gothenburg, Göteborgs universitet, Institutionen för geovetenskaper, 2005. Statistical downscaling of climate scenarios over Scandinavia. Climate Research 29, 255268. doi:10.3354/cr029255  Hatfield, J., Takle, G., Grotjahn, R., Holden, P., Izaurralde, R.C., Mader, T., Marshall, E.,  Liverman, D., 2014: Ch. 6: Agriculture. Climate Change Impacts in the United States: The Third National Climate Assessment, J. M. Melillo, Terese (T.C.) Richmond, and G. W. Yohe, Eds., U.S. Global Change Research Program, 150-174. doi:10.7930/J02Z13FR.  Hayhoe, K., VanDorn, J., Croley, T.II, Schlegal, N., Wuebbles, D., 2010. Regional Climate  Change Projections for Chicago and the US Great Lakes. Journal of Great Lakes Research 36, 721. doi:10.1016/j.jglr.2010.03.012  Hoffmann, H., Rath, T., 2013. Future Bloom and Blossom Frost Risk for Malus domestica  Considering Climate Model and Impact Model Uncertainties. PLOS ONE 8, e75033. doi:10.1371/journal.pone.0075033  Hur, J., Ahn, J., 2015. The change of firstflowering date over South Korea projected from  downscaled IPCC AR5 simulation: peach and pear. International Journal of Climatology 35, 19261937. doi:10.1002/joc.4098  Jackson, M.T., 1966. Effects of Microclimate on Spring Flowering Phenology. Ecology 47, 407 415. doi:10.2307/1932980  Kanamitsu, M., Ebisuzaki, W., Woollen, J., Yang, S.K., Hnilo, J.J., Fiorino, M., Potter, G.L.,  138  2002. NCEP-DOE AMIP-II reanalysis (R-2). Bulletin of the American Meteorological Society 83, 16311643. doi:10.1175/BAMS-83-11-1631  Kaukoranta, T., Tahvonen, R., Ylamaki, A., 2010. Climatic potential and risks for apple growing  by 2040. Agricultural and Food Science 19, 144159. doi:10.2137/145960610791542352  Kunkel, K.E, Stevens, L.E., Stevens, S.E., Sun, L., Janssen, E., Wuebbles, D., Hilberg, S.D.,  Timlin, M.S., Stoecker, L., Westcott, N.E., Dobson, J.G., 2013: Regional Climate Trends and Scenarios for the U.S. National Climate Assessment. Part 3. Climate of the Midwest U.S., NOAA Technical Report NESDIS 142-3, 95 pp.   Frost Risk in the Sour Cherry (Prunus Cerasus L.) Blooming Period in 19852010. Bulletin of Geography. Physical Geography Series 6, 99115. doi:10.2478/bgeo-2013-0006  Ladányi, M., Persely, S., Szabó, T., Szabó, Z., Soltész, M., Nyéki, J., 2010. Climatic indicator  analysis of blooming time for sour cherries. Int. J. Horticultural Science 16, 1116.  Legave, J., Blanke, M., Christen, D., Giovannini, D., Mathieu, V., Oger, R., 2013. A  comprehensive overview of the spatial and temporal variability of apple bud dormancy release and blooming phenology in Western Europe. International Journal of Biometeorology 57, 317331. doi:10.1007/s00484-012-0551-9  Lindkvist, L., Gustavsson, T., Bogren, J., 2000. A frost assessment method for mountainous  areas. Agricultural and Forest Meteorology 102, 5167. doi:10.1016/S0168-1923(99)00087-8  Liszewska, M., Konca-Kedzierska, K., Winkler, J., Tan, P.-N., Abraham, Z., Perdinan, Zhong,  S., Smialecka, E., 2012. Analysis of climate change in selected regions with tart cherries orchards in Central Europe and Michigan (USA) based on regional climate model simulations. Presented at the International Geographical Congress, IGC-2012.  Lobell, D., Field, C., 2011. California perennial crops in a changing climate. Climatic Change  109, 317333. doi:10.1007/s10584-011-0303-6  Logan, J., Mueller, M.A., Searcy, M.J., 2000. Microclimates, peach bud phenology, and freeze  risks in a topographically diverse orchard. HortTechnology 10, 337340.  Longstroth, M., 2002. Freeze damage depends on tree fruit stage of development. Michigan State  University Extension. URL http://msue.anr.msu.edu/news/freeze_damage_depends_on_tree_fruit_stage_of_development  Luedeling, E., Zhang, M., Girvetz, E.H., 2009. Climatic Changes Lead to Declining Winter Chill  139  for Fruit and Nut Trees in California during 19502099. PLOS ONE 4, e6166. doi:10.1371/journal.pone.0006166  Madelin, M., Beltrando, G., 2005. Spatial interpolation-based mapping of the spring frost hazard  in the Champagne vineyards. Meteorological Applications 12, 5156. doi:10.1017/S1350482705001568  MDARD, 2014. Michigan Agriculture Facts & Figures. Michigan Department of Agriculture and  Rural Development.  Mearns, L.O., Gutowski, W., Jones, R., Leung, R., McGinnis, S., Nunes, A., Qian, Y., 2009. A  Regional Climate Change Assessment Program for North America. Eos Trans. AGU 90, 311311. doi:10.1029/2009EO360002  Mearns, L.O., McGinnis, S., Arritt, R., Biner, S., Duffy, P., Gutowski, W., Held, I., Jones, R.,  Leung, R., Nunes, A., Snyder, M., Caya, D., Correia, J., Flory, D., Herzmann, D., Laprise, R., 2007. The North American Regional Climate Change Assessment Program dataset. National Center for Atmospheric Research Earth System Grid data portal, Boulder, CO. doi:10.5065/D6RN35ST  Mearns, L.O., Sain, S., Leung, L.R., Bukovsky, M.S., McGinnis, S., Biner, S., Caya, D., Arritt,  R.W., Gutowski, W., Takle, E., Snyder, M., Jones, R.G., Nunes, A.M.B., Tucker, S., Herzmann, D., McDaniel, L., Sloan, L., 2013. Climate change projections of the North American Regional Climate Change Assessment Program (NARCCAP). Climatic Change 120, 965975. doi:10.1007/s10584-013-0831-3  Menzel, A., Sparks, T.H., Estrella, N., Koch, E., Aasa, A., Ahas, R., Alm-Kübler, K., Bissolli, P.,  Defila, C., Donnelly, A., Filella, Y., Jatczak, K., Måge, F., Mestre, A., Nordli, Ø., r, H., Striz, M., Susnik, A., Van Vliet, A.J.H., Wielgolaski, F.-E., Zach, S., Zust, A., 2006. European phenological response to climate change matches the warming pattern. Global Change Biology 12, 19691976. doi:10.1111/j.1365-2486.2006.01193.x  Mitchell, K.E., Lohmann, D., Houser, P.R., Wood, E.F., Schaake, J.C., Robock, A., Cosgrove,  B.A., Sheffield, J., Duan, Q., Luo, L., Higgins, R.W., Pinker, R.T., Tarpley, J.D., Lettenmaier, D.P., Marshall, C.H., Entin, J.K., Pan, M., Shi, W., Koren, V., Meng, J., Ramsay, B.H., Bailey, A.A., 2004. The multi-institution North American Land Data Assimilation System (NLDAS): Utilizing multiple GCIP products and partners in a continental distributed hydrological modeling system. Journal of Geophysical Research: Atmospheres 109. doi:10.1029/2003JD003823  Moghadam, E.G., Hosseini, P., Mokhtarian, A., 2009. Blooming phenology and self- incompatibility of some commercial cherry (Prunus avium L.) cultivars in Iran. Scientia horticulturae 123, 2933. doi:10.1016/j.scienta.2009.07.013  140  Molitor, D., Caffarra, A., Sinigoj, P., Pertot, I., Hoffmann, L., Junk, J., 2014. Late frost damage  risk for viticulture under future climate conditions: a case study for the Luxembourgish winegrowing region. Australian Journal of Grape and Wine Research 20, 160168. doi:10.1111/ajgw.12059  Mosedale, J.R., Wilson, R.J., Maclean, I.M.D., 2015. Climate Change and Crop Exposure to  Adverse Weather: Changes to Frost Risk and Grapevine Flowering Conditions. PLOS ONE 10, e0141218. doi:10.1371/journal.pone.0141218  Moss, R.H., Edmonds, J.A., Hibbard, K.A., Manning, M.R., Rose, S.K., van Vuuren, D.P.,  Carter, T.R., Emori, S., Kainuma, M., Kram, T., Meehl, G.A., Mitchell, J.F.B., Nakicenovic, N., Riahi, K., Smith, S.J., Stouffer, R.J., Thomson, A.M., Weyant, J.P., Wilbanks, T.J., 2010. The next generation of scenarios for climate change research and assessment. Nature 463, 747756.  NASS, 2015. Cherry Production. National Agricultural Statistics Service, Agricultural Statistics  Board, U.S. Department of Agriculture.  NASS, 2012. Cherry Production. National Agricultural Statistics Service, Agricultural Statistics  Board, U.S. Department of Agriculture.  NASS, 2002. Cherry Production. National Agricultural Statistics Service, Agricultural Statistics  Board, U.S. Department of Agriculture.  Öztekin, T., 2008. Analysis of frost damage risks in peach orchards around Tokat, Turkey.  International J. Natural and Engineering Sciences 2, 4551.  Pagter, M., Arora, R., 2013. Winter survival and deacclimation of perennials under warming  climate: physiological perspectives. Physiologia plantarum 147, 7587.  Perry, K.B., 1998. Basics of frost and freeze protection for horticultural crops. HortTechnology  8, 1015.  Pettit, A.N., 1979. A non-parametric approach to the change-point detection. Appl. Stat 28, 126 135.  PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, created 4 Feb  2004.  Pryor, S., Barthelmie, R., Schoof, J., 2013. High-resolution projections of climate-related risks  for the Midwestern USA. Climate Research 56, 6179. doi:10.3354/cr01143  Pryor, S. C., Scavia, D., Downer, C., Gaden, M., Iverson, L., Nordstrom, R., Patz, J., Robertson,  G.P., 2014: Ch. 18: Midwest. Climate Change Impacts in the United States: The Third National Climate Assessment, J. M. Melillo, Terese (T.C.) Richmond, and G. W. Yohe, Eds., U.S. Global Change Research Program, 418-440. doi:10.7930/J0J1012N. 141   Rochette, P., Bélanger, G., Castonguay, Y., Bootsma, A., Mongrain, D., 2004. Climate change  and winter damage to fruit trees in eastern Canada. Canadian Journal of Plant Science 84, 11131125. doi:10.4141/P03-177  Rodrigo, J., 2000. Spring frosts in deciduous fruit trees  morphological damage and flower  hardiness. Scientia Horticulturae 85, 155173. doi:10.1016/S0304-4238(99)00150-8  San Martino, L., Manavella, F.A., García, D.A., Salato, G., 2005. Phenology and fruit quality of  nine sweet cherry cultivars in South Patagonia, in: V International Cherry Symposium 795. pp. 841848.  Scheifinger, H., Menzel, A., Koch, E., Peter, C., 2003. Trends of spring time frost events and  phenological dates in Central Europe. Theoretical and Applied Climatology 74, 4151. doi:10.1007/s00704-002-0704-6  Schultze, S.R., Sabbatini, P., Andresen, J.A., 2014. Spatial and Temporal Study of Climatic  Variability on Grape Production in Southwestern Michigan. American Journal of Enology and Viticulture 65, 179188. doi:10.5344/ajev.2013.13063  Snyder, R.L., De Melo-Abreu, P.J., 2005. Frost protection: fundamentals, practice and  economics. Food and Agriculture Organization of the United Nations.  Special Report on Emissions Scenarios, 2000. A Special Report of Working Group III of the  Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge.  Taylor, K.E., Stouffer, R.J., Meehl, G.A., 2011. An Overview of CMIP5 and the Experiment  Design. Bull. Amer. Meteor. Soc. 93, 485498. doi:10.1175/BAMS-D-11-00094.1  Themeßl, M.J., Gobiet, A., Heinrich, G., 2012. Empirical-statistical downscaling and error  correction of regional climate models and its impact on the climate change signal. Climatic Change 112, 449468.  Themeßl, M.J., Gobiet, A., Leuprecht, A., 2011. Empirical-statistical downscaling and error  correction of daily precipitation from regional climate models. Int. J. Climatol. 31, 15301544. doi:10.1002/joc.2168  Thornton, P.K., Ericksen, P.J., Herrero, M., Challinor, A.J., 2014. Climate variability and  vulnerability to climate change: a review. Global Change Biology 20, 33133328. doi:10.1111/gcb.12581  University Corporation for Atmospheric Research (UCAR), 2007. Caveats for Users. North  American Regional Climate Change Assessment Program. URL http://www.narccap.ucar.edu/about/caveats.html (accessed 5.25.16).  Vavrus, S., Walsh, J.E., Chapman, W.L., Portis, D., 2006. The behavior of extreme cold air  142  outbreaks under greenhouse warming. International Journal of Climatology 26, 11331147. doi:10.1002/joc.1301  Vitasse, Y., Lenz, A., Körner, C., 2014. The interaction between freezing tolerance and  phenology in temperate deciduous trees. Front Plant Sci 5. doi:10.3389/fpls.2014.00541  Von Neumann, J., 1941. Distribution of the ratio of the mean square successive difference to the  variance. The Annals of Mathematical Statistics 12, 367395.  Walsh, J., Wuebbles, D., Hayhoe, K., Kossin, J., Kunkel, K., Stephens, G., Thorne, P., Vose, R.,  Wehner, M., Willis, J., Anderson, D., Kharin, V., Knutson, T., Landerer, F., Lenton, T., Kennedy, J., and Somerville, R., 2014: Appendix 3: Climate Science Supplement. Climate Change Impacts in the United States: The Third National Climate Assessment, J. M. Melillo, Terese (T.C.) Richmond, and G. W. Yohe, Eds., U.S. Global Change Research Program, 735-789. doi:10.7930/J0KS6PHH.  Watanabe, S., Kanae, S., Seto, S., Yeh, P.J.F., Hirabayashi, Y., Oki, T., 2012. Intercomparison of  bias-correction methods for monthly temperature and precipitation simulated by multiple climate models. Journal of Geophysical Research. Atmospheres 117. doi:10.1029/2012JD018192  Wijngaard, J.B., Klein Tank, A.M.G., Können, G.P., 2003. Homogeneity of 20th century  European daily temperature and precipitation series. International Journal of Climatology 23, 679692.  Wilby, R.L., Charles, S.P., Zorita, E., Timbal, B., Whetton, P., Mearns, L.O., 2004. Guidelines  for use of climate scenarios developed from statistical downscaling methods.  Winkler, J., 2016. (In press) Embracing complexity and uncertainty. Annals of the Association of  American Geographers.  Winkler, J., Andresen, J.A., Bisanz, J.M., Guentchev, G., Nugent, J., Piromsopa, K., Rothwell,  N., Zavalloni, C., Clark, J., Min, H.K., Pollyea, A., Prawiranata, H., 2012b. Ch. 8 Sour Cherry Industry: Vulnerability to Climate Variability and Change, in: Pryor, S.C. (Ed.), Climate Change in the Midwest: Impacts, Risks, Vulnerability, and Adaptation. Indiana University Press, Bloomington, IN, pp. 6981.  Winkler, J.A., Arritt, R.W., Pryor, S.C. 2012a: Climate Projections for the Midwest: Availability,  Interpretation and Synthesis. In: U.S. National Climate Assessment Midwest Technical Input Report. J. Winkler, J. Andresen, J. Hatfield, D. Bidwell, and D. Brown, coordinators. Available from the Great Lakes Integrated Sciences and Assessment (GLISA) Center, http://glisa.msu.edu/docs/NCA/MTIT_Future.pdf. US National Climate Assessment Midwest Technical Input Report.  Winkler, J.A., Cinderich, A.B., Ddumba, S.D., Doubler, D., Nikolic, J., Perdinan, Pollyea, A.M.,  143  Young, D.R., Zavalloni, C., 2013. 2.04 - Understanding the Impacts of Climate on Perennial Crops, in: Climate Vulnerability: Understanding and Addressing Threats to Essential Resources. Elsevier, pp. 3749.   Winkler, J.A., Guentchev, G.S., Liszewska, M., Perdinan, Tan, P.-N., 2011a. Climate Scenario  Development and Applications for Local/Regional Climate Change Impact Assessments: An Overview for the Non-Climate Scientist: Part II: Considerations When Using Climate Change Scenarios. Geography Compass 5, 301328. doi:10.1111/j.1749-8198.2011.00426.x  Winkler, J.A., Guentchev, G.S., Perdinan, Tan, P.-N., Zhong, S., Liszewska, M., Abraham, Z.,  Local/Regional Climate Change Impact Assessments: An Overview for the Non-Climate Scientist: Part I: Scenario Development Using Downscaling Methods. Geography Compass 5, 275300. doi:10.1111/j.1749-8198.2011.00425.x  Winkler, J.A., Thornsbury, S., Artavia, M., Chmielewski, F.-M., Kirschke, D., Lee, S.,  Liszewska, M., Loveridge, S., Tan, P.-N., Zhong, S., Andresen, J.A., Black, J.R., Kurlus, R., Nizalov, D., Olynk, N., Ustrnul, Z., Zavalloni, C., Bisanz, J.M., Bujdosó, G., Fusina, L., Henniges, Y., Hilsendegen, P., Lar, K., Malarzewski, L., Moeller, T., Murmylo, R., Niedzwiedz, T., Nizalova, O., Prawiranata, H., Rothwell, N., Ravensway, J. van, Witzke, H. von, Woods, M., 2010. A conceptual framework for multi-regional climate change assessments for international market systems with long-term investments. Climatic Change 103, 445470. doi:10.1007/s10584-009-9781-1  Wuebbles, D., Hayhoe, K., 2004. Climate Change Projections for the United States Midwest.  Mitigation and Adaptation Strategies for Global Change 9, 335363. doi:10.1023/B:MITI.0000038843.73424.de  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. Journal of Geophysical Research: Atmospheres 117. doi:10.1029/2011JD016048  Zavalloni, C., Andresen, J.A., Flore, J.A., 2006. Phenological Models of Flower Bud Stages and  -day Accumulation. J. Amer. Soc. Hort. Sci. 131, 601607.  Zhao, M., Peng, C., Xiang, W., Deng, X., Tian, D., Zhou, X., Yu, G., He, H., Zhao, Z., 2013.  Plant phenological modeling and its application in global climate change research: overview and future challenges. Environmental Reviews 21, 114. doi:10.1139/er-2012-0036  144