Warm winters challenge the cultivation of temperate species in South America—a spatial analysis of chill accumulation

Winter chill accumulation plays a crucial role in determining the moment of bud burst in temperate fruit and nut trees, and insufficient chill can greatly limit yield potentials. To assess future cultivation options for such species in South America, we estimated winter chill through a spatial analysis. We used historical data (1980–2017) from 158 weather stations to calibrate a weather generator and produce temperature data for 10 historical and 60 future scenarios. We computed Safe Winter Chill (SWC, corresponding to the 10th quantile of a chill distribution) for the period 1980–2017 and for all historical and future weather scenarios and developed a framework to interpolate SWC for the continent using the Kriging method. To improve the interpolation, we applied a 3D correction model based on two co-variables (means of daily temperature extremes in July). Our results suggest important chill declines in southern Brazil and central Chile. By 2050 under the global warming scenario RCP4.5 (Representative Concentration Pathway), absolute SWC in these regions may reach a median of 18.7 and 39.6 Chill Portions (CP), respectively. Projections are most alarming for a strong global warming scenario (RCP8.5). In southern South America, adequate SWC levels of about 60 CP may be expected even under the RCP8.5 scenario. Our results highlight the need for climate change adaptation measures to secure temperate fruit production in important growing regions of South America. The procedure we developed may help farmers and practitioners across South America estimate future SWC to adapt their orchards to future challenges.


Introduction
Woody species that evolved in temperate climate regions enter a dormancy phase in winter that allows these plants to withstand cold and freezing temperatures. During dormancy, all visible growth is suspended until the onset of favorable conditions in spring (Fadón et al. Previous studies have assessed the impacts of climate change on agroclimatic metrics in different regions and climates. Projections suggest a major chill decline in Mediterranean climate areas such as the Central Valley of California, Tunisia, and some areas of Australia, Chile, Oman, and Spain (Benmoussa et al. 2020;Buerkert et al. 2020;Darbyshire et al. 2016;Fernandez et al. 2020b;Rodríguez et al. 2019). In other regions such as Patagonia in Argentina, projections are less alarming regarding chill decline and suggest increasing levels of heat during the growing season, which may improve conditions for fruit production (del Barrio et al. 2021). Many of the studies mentioned above assessed the impacts of climate change on specific locations by using weather data recorded in situ. Although this approach produces relevant information on historical and future climatic trends, analyzing chill levels in specific locations may restrict the applicability of possible recommendations to only the sites where data were recorded, with farmers outside these locations remaining without adequate support. Recent studies by Benmoussa et al. (2020) and Rodríguez et al. (2019) have addressed this limitation by using spatial interpolation of site-specific chill estimates. In addition to spatial interpolation of simulation data, Benmoussa et al. (2020) accounted for the elevation profile of the region to generate a chill surface for Tunisia. The elevation profile in a region has previously been demonstrated to be an important factor modulating chill accumulation in mountainous regions (Alburquerque et al. 2008;Buerkert et al. 2020). Besides the elevation profile, other proxies of chill accumulation that are related to specific geographic or climatic characteristics of a region may help interpolate chill levels for large areas where the scarcity of weather data recorded in situ is a limiting factor.
South American agriculture contributes to local and global food security. Some regions of Peru, Brazil, Argentina, and Chile are characterized by substantial production of temperate fruits, which places these countries among the most relevant exporters of fresh fruit at the global scale (CAFI 2019;ODEPA 2017). To secure the cultivation of temperate species among the different climates observed on the continent, decision makers, farmers, and orchard managers must be provided with adequate recommendations based on up-todate knowledge. At present, the only available map of chill accumulation that includes the whole region was produced by  in a global analysis implemented about a decade ago. In their study,  used the A2, A1B, and B1 scenarios of greenhouse gas emissions for temperature projections (with A1B and B1 being comparable to RCP6.0 and RCP4.5, respectively according to the US Global Change Research Program (Global Change 2021)). A regional analysis for southern Brazil using the Chilling Hours model (Weinberger 1950) was reported by Wrege et al. (2010). While the results of the study by  and Wrege et al. (2010) provided useful insights at a global and regional scale, respectively, a continental analysis with updated warming scenarios is likely to improve our understanding of the potential risks regarding chill levels for farms cultivating temperate trees in different regions of South America.
The main aim of this work was to estimate and map historical and future chill levels for South America. We quantified chill accumulation in terms of Chill Portions (CP) according to the Dynamic Model (Erez et al. 1990;Fishman et al. 1987aFishman et al. , 1987b. Among other commonly used chill modeling approaches (e.g., the Chilling Hours model cited by Weinberger (1950) or the Utah Model developed by Richardson et al. (1974)), the Dynamic Model has emerged as the most reliable model for chill estimation in mild-winter climates due to its biologically credible structure (Fernandez et al. 2020c;Zhang and Taylor 2011). In brief, the Dynamic Model states that chill accumulates in a two-step process that is modulated by cold and subsequent moderate temperatures (Luedeling 2012). Using this chill model, we estimated chill levels for 158 locations across South America using data recorded in situ between 1980 and 2017 by weather stations across the continent.
For the past, we computed chill based on actual observations as well as based on historical weather scenarios (scenarios that represent typical climatic conditions for specific past periods). To obtain future chill levels, we used temperature projections by an ensemble of fifteen General Circulation Models (GCMs), for two RCP scenarios (RCP4.5 and RCP8.5), and for two time horizons (2050 and 2085). We then interpolated chill levels for each scenario for the whole continent using the Kriging method. We produced the final maps by correcting the linear interpolation using a model that considers the relationship between average daily minimum and maximum temperature in July and observed chill, as recorded at the weather stations.

Selection of weather stations and data collection
We selected 158 weather stations located across South America as sources of temperature data ( Fig. 1; Table S1). To select appropriate stations, we first listed 1500 weather  . We filtered the initial set of weather stations to select those with data between January 01, 1980 and December 31, 2017. This resulted in 332 weather stations, for which we downloaded weather records from the GSOD database. After downloading, we determined the records' coverage of minimum and maximum temperatures for the intended period of this study (i.e., the chilling season between May 1 and August 31) and selected all stations whose records were at least 90% complete for at least one of the two metrics. This reduced the number of weather stations to 105. In order to increase the reliability of the analysis, we included weather data from different sources. We accessed a Chilean network maintained by the Center for Climate and Resilience Research ([CR] 2 ) through functions contained in the dormancyR package (Fernandez 2020). Out of 909 available weather stations across the country, we selected 58 stations that had at least 90% data coverage for the winter period of 1980-2017. Additionally, we obtained data for 5 weather stations located in Patagonia, Argentina that fulfilled the data completeness requirement. These data were collected by the National Weather Service (SMN, https:// www. smn. gob. ar) as well as the National Institute of Agricultural Technology (INTA, https:// www. argen tina. gob. ar/ inta) and had been used in an earlier study by del Barrio et al. (2021).

Data preparation and gap filling
The data downloaded from the databases showed a maximum of 9.9% of missing values for daily minimum temperatures and 16.5% of missing values for daily maximum temperatures. On average, the percentage of missing values among weather stations was ~ 3% (about 145 days out of 4674) for both minimum and maximum temperatures. To fill the gaps with auxiliary data sources, we applied the procedure described in Fernandez et al. (2020c), with some modifications. In brief, we identified the 40 closest weather stations to each primary weather station that showed at least one gap. For each weather station in this auxiliary list of 40 data sources, we computed the difference relative to the primary location for all days that had records in both datasets. This procedure was implemented separately for minimum and maximum temperatures and for each month of the year. We then computed the monthly mean difference as well as the monthly standard deviation for the resulting values. We used the mean difference to adjust the data in the auxiliary dataset before filling the gaps in the primary dataset. To avoid including non-representative data, we set a threshold of 4 °C for both the monthly mean difference and standard deviation to filter the auxiliary weather stations. We applied this procedure until filling all gaps or until no more auxiliary weather stations were available.
The procedure described above resulted in complete records for most weather stations. Nonetheless, a relevant share (75 weather stations) still showed missing days after applying the patching procedure. In these cases, we filled the remaining gaps through linear interpolation (Luedeling 2018). On average, ~ 30 days were linearly interpolated in each weather station (ranging from 1 to 214 days). For most of the stations (52 locations), the interpolated days represented periods outside the chill accumulation season, which were not relevant for further analyses. However, in 23 out of 75 weather stations, the linear interpolation procedure covered the period we later used for estimating chill accumulation (from May to August). For these stations, only between 1 and 4 days were interpolated, suggesting that our analysis was not greatly influenced by applying linear interpolation to fill data gaps (Luedeling 2018). After applying the gap-filling procedure, ten weather stations were removed from the analysis since these stations showed missing values for more than a month. The steps described above helped us obtain complete records of daily temperature extremes (minimum and maximum records) for 158 locations (supplementary materials Table S1) across South America for a period of 38 years.

Historical and future scenarios
We computed winter chill for all years on record for each weather station. Computing chill accumulation from observed weather, however, may obscure long-term trends, which are crucial for assessing future climate risks. Each computation of chill for a given year is essentially the result of a random draw from an unknown distribution of plausible temperature outcomes. This logic underlies the functioning of weather generators, which derive these unknown distributions from observed weather records and offer capabilities to draw such random samples. To characterize agroclimatic conditions in the past, we therefore generated historical weather scenarios based on the data recorded on-site. As reported in previous studies (Benmoussa et al. 2020;Buerkert et al. 2020), we used the RMAW-GEN weather generator (Cordano and Eccel 2019) to produce 100 replicates of synthetic weather and characterize the typical distributions of weather conditions for particular years in the past. Weather generators are statistical models that simulate realistic random sequences of atmospheric variables such as temperature, precipitation, and wind (Wilks and Wilby 1999) based on a calibration dataset. As inputs for this tool, we used monthly means of daily minimum and maximum temperatures for each year of interest, computed via a 15-year running mean function applied to data for all years on record. This allowed us to closely follow patterns observed in the historical data. We used this procedure to generate scenarios for 1981, 1985, 1989, 1993, 1997, 2001, 2005, 2009, 2013, and 2017. For each of these years, the outputs produced by the weather generator represented plausible distributions of temperature conditions that are characteristic of the respective year.
We used the same procedure described above to generate future climate scenarios. As input for the weather generator in this case, we used monthly means of daily temperature extremes downloaded from the ClimateWizard ) database maintained by the International Center for Tropical Agriculture (CIAT) using an application programming interface (https:// github. com/ CIAT-DAPA/ clima te_ wizard_ api). ClimateWizard contains temperature projections generated by 15 General Circulation Models (GCMs, supplementary materials Table S2) for two Representative Concentration Pathway scenarios (RCP4.5 and RCP8.5) and for a number of time horizons. For this study, we selected the periods 2035-2065 and 2070-2100, denoted by their central years 2050 and 2085. Since the ClimateWizard database contains temperature projections relative to the period 1950-2005 (by default), we implemented a baseline adjustment procedure to make the estimations comparable to our historical period. We first generated two absolute temperature scenarios based on the weather data recorded on-site. These temperature scenarios represented the conditions (mean minimum and maximum temperatures for each month) of the central years of the periods 1980-2017 (historical weather data) and 1980-2005 (baseline selected in the ClimateWizard database), respectively. We then obtained a relative temperature scenario by computing the monthly difference for each metric between these two absolute scenarios. We used the resulting relative scenario to adjust the temperature scenario obtained from ClimateWizard. The final output was used to feed the weather generator and produce 100 years of weather records based on projections by 15 GCMs.

Winter chill estimation
We computed winter chill for (i) all years on record, (ii) all synthetic years for each historical scenario, and (iii) all synthetic years for each combination of RCP, time horizon, and GCM. For the sake of standardization and comparison among past and future scenarios as well as with results from previous studies, we selected fixed dates for the start and end of the chill accumulation season. Following dates used in previous studies (see Fernandez et al. 2020c for example), we defined the chill accumulation season as between May 1 and August 31 for each year. We used the Dynamic Model (Erez et al. 1990;Fishman et al. 1987aFishman et al. , 1987b to compute chill, since this model has emerged from a number of model comparison studies as the most credible tool for explaining the process of chill accumulation in temperate trees (Fernandez et al. 2020c;Zhang and Taylor 2011). Because this model requires hourly temperatures as input, we derived hourly records from daily temperature extremes and the latitude of each weather station. To this end, we estimated temperatures based on an idealized daily temperature curve consisting of a sine function for daytime warming and a logarithmic decay function for nighttime cooling (Almorox et al. 2005;Linvill 1990).
We summarized the chill distributions for each weather station and scenario according to the Safe Winter Chill concept (SWC, Luedeling et al. 2009). Safe Winter Chill is the 10th percentile of the chill distribution, which represents the chill level that is likely to be met in 90% of years (Luedeling et al. 2009). Regarding future chill levels computed across the 15 GCMs, we summarized the Safe Winter Chill values according to 'optimistic', 'intermediate', and 'pessimistic' GCM projections. These climate model classes represented the 85th, 50th, and 15th quantile of the chill distribution across all GCMs (15 values) for each combination of RCP and time horizon, respectively.

Predictors of Safe Winter Chill
To estimate winter chill accumulation in areas lacking observational data, we identified geographic and environmental variables correlated with SWC. We first assessed the linear relationship between SWC observed on site and the variables of interest. To estimate the potential effect of latitude on such a relationship, we divided the list of weather stations into 12 groups from 5° N to 55° S. We tested the distance to the ocean (in km), the elevation above sea level (in m), and the average of daily minimum and maximum temperatures in July (in °C). Based on the results of this analysis, we decided to use the average of daily minimum and maximum temperatures in July as predictors of SWC for the spatial interpolation.

Spatial chill interpolation
We developed a methodological framework (Fig. 2) for spatial interpolation of SWC, building on methods described by Benmoussa et al. (2020). We first interpolated an SWC surface using an ordinary Kriging approach. Since such Kriging surfaces are only based on data for the particular locations used for the interpolation, they may contain large errors for many other regions. To account for the considerable variation in climate and geography across the continent, which is known to affect chill accumulation, we established a correction procedure based on two co-variables with known spatial distributions. Our correction variables, which may serve as proxies for chill accumulation, were the means of daily minimum and maximum temperatures in July. We downloaded historical gridded data  for these co-variables with a spatial resolution of 30 arcseconds (~ 1 km 2 ) from the WorldClim database (Fick and Hijmans 2017). In addition to these gridded data, we used ordinary Kriging to produce interpolated surfaces for the same variables recorded in situ by the weather stations, which were used to characterize the temperature conditions that were represented by the original Kriged SWC surface. To include only representative data, we excluded weather stations for which the difference between proxy variable values from WorldClim and the weather stations was greater than 2 °C (see supplementary materials Fig. S1).
For each scenario, we created a 3D correction model by generating a surface contour plot relating SWC to mean daily minimum and maximum temperatures in July computed for the weather stations (Fig. 2). To create the 3D correction model, we assumed that observed SWC values are a function of the daily minimum and maximum temperatures in July and implemented ordinary Kriging to predict a continuous surface. We implemented this procedure using functions within the fields package for R (Nychka et al. 2017). We interpolated a similar surface relating SWC to values of the co-variables for the same locations from the WorldClim database. For each point on the map, we then subtracted the SWC value from the contour plot based on weather station data from the contour plot based on WorldClim data. The result was a correction SWC map that we added to the original SWC values obtained through ordinary Kriging (Fig. 2). We note that in this correction Fig. 2 Framework for spatial interpolation of winter chill across South America. In the first stage, we applied ordinary Kriging to data for weather station locations to obtain surfaces of Safe Winter Chill (SWC) and mean daily minimum and maximum temperatures in July using data recorded in situ ('Chill map' and 'Correction variable maps' nodes in the figure). In a second step, we obtained gridded data for our correction co-variables from the WorldClim database ('Correction variable maps' node). Using data from the weather stations, we then generated a '3D correction model' based on the relationship between monthly means of daily temperature extremes (minimum and maximum records) and observed SWC. We used this model to estimate SWC values from the co-variables (mean daily minimum and maximum temperatures in July) from both data sources (WorldClim and weather stations) to obtain the SWC correction map ('Correction map' node). Finally, we applied this SWC correction map to the original Kriged map to produce the final SWC maps for all scenarios 1 3 procedure, we assumed that the error in the original SWC Kriging is related to the difference between co-variables from both databases (WorldClim and on-site). We set the correction procedure to generate chill values ≥ 0 considering that the Dynamic Model returns positive chill values only. Additionally, we excluded all regions (using a hatching layer to mask them in all maps) for which the combination of mean daily minimum and maximum temperatures in July had no observed (on-site) values of SWC in the 3D correction model (see supplementary materials Fig. S2).

Cross-validation
We implemented a cross-validation procedure to assess the reliability of our spatial interpolation approach. After excluding all weather stations that featured a difference > 2 °C compared with the WorldClim data for mean daily minimum or maximum temperatures in July, we randomly split the set of weather stations into 7 groups. Each of the seven groups was once treated as a validation data set, while the remaining data sets were combined into a training set to compile the final SWC maps following our interpolation framework (Fig. 2). We then compared the values for SWC computed with the training data set to the values observed in the validation set to obtain the residual values (observed minus predicted). Since the grouping can affect the calculation of residuals, we repeated the crossvalidation procedure 5 times, and computed the mean residual across these 5 repetitions. We iterated this process over the observed and simulated historical scenarios as well as the 'pessimistic', 'intermediate', and 'optimistic' climate model classes by 2050 and 2085 under RCP4.5 and RCP8.5. We summarized the results of the cross-validation procedure by computing the median and standard deviation of residuals across scenarios.

Data processing, analyses, and figure preparation
All data processing, analyses, and figure preparation were implemented in the R programming environment (R Core Team 2020). Analyses related to agroclimatic metrics were performed with functions of the chillR package . For figure preparation, we used the ggplot2 (Wickham 2016) and tmap (Tennekes 2018) packages. Scripts for reproduction of all data processing steps (including figure generation) are available in a public GitHub repository (https:// github. com/ Eduar doFer nande zC/ chill_ south_ ameri ca).

Results
To facilitate visualization of results, we only included figures representing key scenarios (e.g., RCP4.5-2050 and RCP8.5-2085 for future conditions) in the main text. For all combinations of RCP scenario, time horizon, and climate model class, see our supplementary materials.

Predictors of SWC
The results of our analysis suggest that the means of daily minimum and maximum temperatures in July may serve as the best predictors of winter chill accumulation in South America when compared with elevation and distance from the ocean (Fig. S3). When assuming a linear relationship between SWC and these variables, the overall coefficient of determination of the regressions showed values of 0.674 for the mean minimum and 0.796 for the mean maximum temperature. We detected much weaker correlations between SWC and distance to the ocean (R 2 of 0.169) and elevation above sea level (R 2 of 0.512). When including the interaction between the variables and the latitude group, R 2 values tended to increase considerably. For the means of daily minimum and maximum temperature, this correlation metric increased to 0.905 and 0.933, respectively. For the remaining variables (distance to the ocean and elevation), the R 2 increased to 0.857 and 0.880.

Cross-validation
Our cross-validation approach resulted in rather small residuals across weather stations, but also identified some notable outliers (Fig. 3). Overall, the median residual was −0.85 CP, considering all weather stations showing residual values after the cross-validation procedure. It should be noted that the procedure did not produce residual values for 17 weather stations. This is explained by these weather stations falling outside the boundaries of the 3D correction model (created with the training dataset) when included in the evaluation dataset. Median residuals ranged from −3.86 CP (Q25%) to 0.93 CP (Q75%) and showed maximum values of −32.1 CP (Teniente Alejandro Velasco Astete International; 13.5° S; 71.9° W; 3310.1 m.a.s.l.) for the lower bound and 7.16 CP (Digua Embalse; 36.3° S; 71.6° W; 390.0 m.a.s.l.) for the upper bound. The standard deviation of the residuals ranged from 1.53 CP (Q25%) to 3.70 CP (Q75%) and showed a median of 2.54 CP across weather stations. The maximum standard deviation was 11.6 CP, identified for the weather station Teniente Alejandro Velasco Astete International (13.5° S; 71.9° W) located at 3310.1 m.a.s.l..
Comparison of SWC calculations for the simulated historical scenarios for 1981 and 2017 indicated that, across the whole continent, changes ranged from −34.3 CP to 16.5 CP (Fig. 4). Although the absolute change was relatively small for most locations, we estimated considerable chill declines near the Andes in Colombia (7-2° N; 72-77° W), in the Central Valley of Chile, and in northern Patagonia, Argentina. In these regions, average SWC change reached −7.9 CP (−14.8 CP for Q5% and 0 CP for Q95%), −5.7 CP (−8.4 CP for Q5% and −3.4 CP for Q95%), and −4.5 CP (−6.8 CP for Q5% and −1.8 CP for Q95%), respectively. In southern Brazil, the overall change in SWC between

Absolute Safe Winter Chill under future conditions
Our results show comparable SWC values among the climate model classes 'pessimistic' and 'optimistic' (Fig. S4 to S6). For the RCP4.5-2050 scenario, maximum SWC values across South America were about 100 CP for the 'pessimistic' class of climate projections, and slightly lower (99.4 CP) for the 'optimistic' class. Conversely, slightly more chill was projected for the 'optimistic' class when compared with the climate model class 'pessimistic' (101.1 CP versus 101.6 CP) under the global warming scenario RCP8.5 by 2085 (Fig.  S4). A considerable SWC decline was observed when comparing the same climate model class (e.g., 'pessimistic') under different RCP scenarios and time horizons (upper versus bottom panels in Fig. S4).
Regarding specific regions of the continent, absolute SWC values are likely to range from 4.9 CP (Q5%) to 28.7 CP (Q95%) and reach a median of 17.0 CP under the 'RCP4.5-2050 pessimistic' scenario in southern Brazil, a major production region of temperate fruit (Figs. S7 to S9). Growers around this region may expect SWC values of 9.9 CP (likely range between 1.4 CP and 17.7 CP-Q5% and Q95%, respectively) under the 'RCP8.5-2085 pessimistic' scenario. Compared with the 'pessimistic' climate model class, projections for an 'optimistic' climate projection show slightly greater SWC values by about 3.0 CP and 5.3 CP for the median for the RCP4.5-2050 and RCP8.5-2085 global warming scenarios, respectively (Fig. S7).
In northern Patagonia, Argentina, SWC is likely to reach a median of 68.9 CP under RCP4.5 by 2050 for the 'pessimistic' climate model class (Fig. S7). For the 'optimistic' class, SWC in this region may increase to 71.3 CP for the same combination of RCP and time horizon. For RCP8.5 by 2085, farmers can expect between 47.1 CP (Q5%) and 82.7 CP (Q95%) according to the 'pessimistic' climate model class and between 56.7 CP (Q5%) and 85.8 CP (Q95%) for the 'optimistic' class.
For the Central Valley of Chile, our projections suggest an SWC range between 25.8 CP (Q5%) and 40.8 CP (95%) for the 'pessimistic' climate model class of the RCP4.5-2050 scenario (median of 36.7 CP). By 2085 and for RCP8.5, chill levels are likely to range between 10.0 CP (Q5%) and 32.4 CP (Q95%), using the same climate model classification (median of 21.2 CP). In the 'optimistic' climate model class, SWC might reach a median of 41.9 CP for the RCP4.5-2050 scenario and 33.3 CP for the RCP8.5-2085 scenario.

Expected changes in chill accumulation compared with the historical period
Compared with the historical scenarios, our results suggest a notable chill decline even for the moderate global warming scenario RCP4.5 (Figs. 5 and S10 to S11). Depending on the region evaluated, chill change ranged between −38.6 CP and 9.9 CP for the 'pessimistic' climate model class of the RCP4.5-2050 scenario and between −29.0 CP and 9.7 CP for the 'optimistic' class for the same RCP and time horizon combination. A much greater SWC decline by about 55 CP and 45 CP might be expected under the RCP8.5-2085 scenario for the 'pessimistic' and 'optimistic' climate model classes, respectively (Fig. 5).
The Central Valley of Chile is likely to emerge as one of the most affected regions across South America regarding winter chill decrease (Fig. 5)  We computed the change using as baseline the median SWC across ten simulated historical scenarios between 1981 and 2017. In the first row ('RCP4.5-2050 pessimistic' and 'RCP4.5-2050 optimistic' panels), we show the expected change for a moderate warming scenario for the near future using two climate model classes. In the second row ('RCP8.5-2085 pessimistic' and'RCP8.5-2085 optimistic' panels), we show SWC changes estimated for a strong warming scenario by the far future using two climate model classes. 'Pessimistic' and 'optimistic' labels represent the 15th and 85th quantiles of the SWC distribution across the 15 General Circulation Models (see supplementary Table S2) we used in this analysis. The grey hatching ('Excluded' in the legend) in the maps indicates regions where the combination of mean daily minimum and maximum temperatures in July fell outside the boundaries of the 3D correction model for the respective scenario. The red crosses represent the location of the weather stations we used in the analysis scenario RCP4.5 for the 'pessimistic' climate model class. In this scenario, the likely decline may range between 16.1 CP (Q5%) and 6.2 CP (Q95%). The 'optimistic' climate model class of the RCP4.5-2050 scenario suggests an SWC decline between 9.9 CP (Q5%) and 3.6 CP (Q95%), with a median of 6.9 CP (Fig. 5). In other regions such as northern Patagonia, Argentina, and southern Brazil, comparable changes are only projected for the RCP8.5-2085 scenario for the 'pessimistic' climate model class.

Discussion
We developed a method for spatial interpolation of Safe Winter Chill on a continental scale. We used mean daily minimum and maximum temperatures in July as co-variables for estimating winter chill accumulation across South America. In a previous study conducted in the Mediterranean climate of Tunisia, Benmoussa et al. (2020) successfully used elevation above sea level as a chill proxy in the interpolation procedure. Similarly, Buerkert et al. (2020) and Alburquerque et al. (2008) highlighted the importance of the elevation of the region in determining chill levels during winter in northern Oman and southeastern Spain. In our study however, elevation above sea level was not a good proxy for winter chill accumulation, probably because of the large extent of the region as well as major intra-regional differences in climatic conditions. In general, the correlation between elevation and winter chill was high at low latitudes (between 5° N and 25° S) but very low at locations south of 30° S (see Fig. S3). Along the same lines, distance to the ocean did not constitute a good proxy for winter chill accumulation across the whole continent. This and other geographic and climatic variables may serve as proxies of chill accumulation in smaller regions or for areas with low variation in topography and climatic conditions. In a highly diverse region (climatically and topographically) such as South America, only monthly means of daily temperature extremes appeared to capture the overall intra-regional variation and exhibit sufficient correlation with SWC. Using these two co-variables to supplement estimates based on on-site records proved effective in obtaining reasonable estimations of winter chill accumulation in South America and holds promise for similar regions where high-resolution datasets are scarce.
We observed large residuals for some weather stations when cross-validating the chill interpolation procedure. Among these, we highlight the weather stations Vilcuya (at 1100 m.a.s.l.) in Chile, Teniente Alejandro Velasco Astete International in Peru (at 3310 m.a.s.l.), and El Dorado International in Colombia (at 2549 m.a.s.l.). The median residuals among repetitions and scenarios for these stations reached −27.6 CP, −32.1 CP, and −20.8 CP, respectively. Large residuals in high-elevation areas may partially be explained by our 3D correction model failing to capture the relationship between the covariables (mean daily temperature extremes in July) and SWC in these regions. Since the 3D correction model is based on on-site data, we believe that more information from these extreme environments may improve the interpolation procedure. Large residuals may also result from particular microclimatic settings at the respective weather stations. Among the three stations with the highest residuals, two are located at high-elevation airports in relatively dense urban areas, whereas the third station is located in a narrow mountain valley. It seems plausible that correction models built without these stations (by them being part of the evaluation dataset) may not be able to represent the special climatic conditions that result from such settings.

3
Our results suggest considerable chill reductions in many parts of the continent even for the moderate global warming scenario. These findings mirror the results of earlier analyses conducted for a number of regions around the globe (Darbyshire et al. 2016;Fernandez et al. 2020b;Rodríguez et al. 2019). The Central Valley of Chile emerges as a highly sensitive region regarding the impacts of global warming on winter chill accumulation. In agreement with previous studies in Chile (Fernandez et al. 2020b(Fernandez et al. , 2020c, orchards located between 30° S and 37° S may experience a reduction in SWC between 6 and 16 CP by 2050 for the RCP4.5 scenario ('pessimistic' climate model class). Analyzing the same combination of RCP and time horizon, absolute SWC levels in this region may reach a median of 36 CP, which may represent a hazardous situation for farmers cultivating high-chill species such as apple and sweet cherry. Whereas apple trees may require between 50 and 90 CP (depending on the cultivar) to successfully overcome dormancy (Delgado et al. 2021;Fernandez et al. 2020a;Parkes et al. 2020), the chill requirement of different sweet cherry cultivars has been reported as between 30 and 86 CP (reviewed in Fadón et al. 2020b). For southern Brazil, Wrege et al. (2010) forecasted absolute chill levels between 300 and 600 Chilling Hours (between 25 and 50 CP using conversion factors across models from Luedeling and Brown 2011) for a scenario considering an increase of 1 °C relative to the period 1961-2000. By 2050 under the RCP4.5 scenario (which projects a likely temperature increase between 1.1 and 2.6 °C), our results suggest a median absolute SWC value of 17.0 CP (28.7 CP for Q95%). Under this scenario, farmers cultivating temperate species in southern Brazil may face severe challenges during the dormancy phase. Providing these farmers with an adequate portfolio of strategies to adapt their orchards to the negative impacts of climate change may become critical for sustaining temperate fruit production.
The application of dormancy-breaking agents is a horticultural practice that may help farmers overcome dormancy-related problems caused by insufficient winter chill. Hydrogen cyanamide, gibberellic acid, and potassium nitrate have successfully been used to promote dormancy release and bud burst in grapevines, Japanese apricots, and pistachios (Perez and Lira 2005;Rahemi and Asghari 2004;Zhuang et al. 2013). For this strategy to be effective, however, a minimum amount of chill (e.g., about 70% of the trees' chill requirement) must have been accumulated by the trees before the chemical is applied (Erez 2000). This basic need limits the application of dormancy-breaking agents to regions with at least some winter chill accumulation, such as central and southern-central Chile and northern Patagonia. In warmer areas (e.g., southern Brazil), farmers may have to develop medium-and long-term strategies to secure temperate fruit production in the future. In this regard, shifting the array of species to low-chill genotypes may offer farmers an opportunity to adapt their orchards to the impacts of climate change. Nonetheless, adopting radical adaptation strategies (e.g., shifting species) is often a decision that farmers hesitate to make because of cultural, market-related, and personal preferences (Nguyen and Drakou 2021).
Developing new low-chill cultivars appears as a promising strategy to help farmers face future challenges while retaining their current production practices for a given species. Breeding new cultivars is, however, a challenging task that requires considering a number of variables such as productivity, consumer preferences, and susceptibility to diseases and growth disorders (Sherman and Beckman 2003). Similarly, breeding programs face the challenge of estimating and reporting the cold requirement of trees, since all metrics to make such estimates were developed more than three decades ago, and none have been fully convincing. Currently available chill models have shown low transferability of estimates among regions (Luedeling 2012;, making reported cultivar-specific chill requirements uninformative for most regions that are different from those where the initial estimation was conducted. In the case of the Dynamic Model, for instance, most farmers and researchers still use the model with the original set of parameters for peach trees reported 30 years ago (Erez et al. 1990;Fishman et al. 1987aFishman et al. , 1987b, although current knowledge suggests that model parameters should be adjusted for each species, and possibly each cultivar (Campoy et al. 2011;Fadón et al. 2020a). Recent studies suggest a need to re-parameterize the Dynamic Model depending on the species and cultivar evaluated (Egea et al. 2020). Similarly,  recently proposed a modelling framework that considers the climatic needs of trees during both the endo-and eco-dormancy phases to predict bloom dates. Including a cultivar-specific re-parameterization, as well as a framework considering the whole dormant phase, is expected to help breeders identify genotypes that are well-adapted to future climatic challenges. Farmers growing cultivars that remain economically viable under future scenarios should be able to secure the viability of temperate fruit production in South America and possibly other warm-winter regions.

Conclusions
Our results raise concerns regarding future chill levels in some areas of South America such as northern Patagonia, southern Brazil, and northern-central Chile. To secure temperate fruit production in these regions, farmers cultivating high-chill species might need to develop new medium-and long-term strategies. Among these, we suggest the cultivation of genotypes with lower chill requirements that are aligned with local conditions under future scenarios. Some regions are likely to be relatively unaffected by chill losses. In southern regions of Chile and Argentina, future chill levels may remain sufficient for overcoming dormancy in temperate fruit species. However, some horticultural practices such as the application of dormancy-breaking agents may play a role in sustaining adequate yields after warm winters. The method we used for interpolating Safe Winter Chill helped estimate chill levels in regions where long-term high-resolution data are scarce. Despite some limitations of the method (low accuracy in high-elevation regions), which indicate a need for additional research, our approach can help provide farmers with reasonable estimates of future chill levels. Estimations of future chill levels help farmers, orchard managers, and relevant stakeholders make decisions on future developments and planning regarding the temperate fruit industry.