A comparison of interpolation methods to predict chill accumulation in a Mediterranean stone fruit production area ( Region de Murcia, SE Spain)

A comparison of interpolation methods to predict chill accumulation in a Mediterranean stone fruit production area (Regi´ on de Murcia, SE Spain) In´ es G´ omez-Ramos a , Manuel Caro b , Juan A. L´ opez b , David Ruiz c , Jose A. Egea c , * a Department of Plant Biology and Ecology, University of Seville, C/ Profesor García Gonz´ alez, s/n, Seville 41012, Spain b Murcia Institute of Agri-Food Research and Development, C/ Mayor, s/n, 30150 Alberca Las Torres, Murcia 30150, Spain c Fruit Breeding Group, Department of Plant Breeding, CEBAS-CSIC, Campus Universitario de Espinardo 25, Murcia 30100, Spain ARTICLE INFO Keywords: Winter chill Spatial interpolation Temperate fruits Cross-validation Stone fruit Climate change ABSTRACT One of the most important agroclimatic variables for stone fruit production is winter chill accumulation. To estimate chill accumulation in locations where climatic data is not recorded, spatial interpolation is necessary. In this study, we compare different interpolation methods for mean and Safe Winter Chill (SWC) in a Mediterranean stone fruit production area (Region of Murcia, SE Spain) using data from 49 climatic stations. To choose the most accurate interpolation method, as its choice may substantially influence the prediction accuracy, the predictive capability of several interpolation methods with different parameterizations (for a total number of 32 instances) was compared through out-of-bag bootstrap cross-validation, concluding that the best ones were Radial Basis Functions applied on the altitude-dependent regression residuals for mean winter chill and the altitude+latitude linear regression for SWC. The incorporation of altitude in the interpolation increased greatly the accuracy of the estimation. In fact, most of the chill accumulation spatial dependency was explained through altitude. The accuracy of the interpolation was not homogeneous across the study area. Chill accumulation in warmer coastal localities was overestimated by all the methods, possibly due to the proximity to the sea, highlighting the importance of microclimatic variables at higher-resolution spatial interpolations. Differences between methods were more notable in higher locations, where distance-only based methods underestimated chill accumulation and methods that consider altitude slightly overestimated it. This study demonstrates the importance of comparing the performance of multiple spatial interpolation methods before applying any for chill accumulation data. 1. Introduction For stone fruit trees the dependence on climatic conditions, and especially on temperature, is essential for their proper development. These trees establish dormancy in late autumn to avoid frost damage in winter or early spring (Fad´ on et al., 2020; Faust et al., 1997; Rohde and Bhalerao, 2007). Two main sub-phases are defined once dormancy has been set: endo-dormancy, where growth does not occur even under suitable environmental conditions, and eco-dormancy where exposure to warm temperatures restarts growth to bloom eventually (Lang et al., 1987). The mechanism to overcome endo-dormancy is based on the concept of chill requirements, a genetic trait that consists in a metric related to cold temperatures within a certain range (Luedeling, 2012). Endo-dormancy is overcome once a certain amount of winter chill is accumulated, fulfilling their chill requirements (Luedeling and Brown, 2011) to minimize the risk of frost. Some authors state that both phases, endo-dormancy and eco-dormancy, can partially overlap and chill and heat accumulation can occur in parallel (Harrington et al., 2010; Pope et al., 2014) or even compensate each other (Kaufmann and Blanke, 2019). Numerous efforts have been made in order to model chill accumulation, due to its importance on yield production. While the most popular chill accumulation models are the so-called Chilling Hours Model (Weinberger, 1950) and the Utah model (Richardson et al., 1974), the most accurate model for temperate fruits in warmer climates is the Dynamic model (Fishman et al., 1987a, 1987b), which is the one recommended by some authors, especially for areas with warm winters (Egea et al., 2021; Fernandez et al., 2020b; Luedeling, 2012). Many of the studies which compare chilling accumulation models, including the Dynamic model, found that its performance was better or as accurate as * Corresponding author. E-mail address: [email protected] (J.A. Egea). Contents lists available at ScienceDirect European Journal of Agronomy journal homepage: www.elsevier.com/locate/eja https://doi.org/10.1016/j.eja.2024.127316 Received 14 March 2024; Received in revised form 7 August 2024; Accepted 7 August 2024 European Journal of Agronomy 160 (2024) 127316 Available online 12 August 2024 1161-0301/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). the rest of the models (Luedeling and Brown, 2011), with some of the studies performed in south-eastern Spain (Alburquerque et al., 2008; Ruiz et al., 2007). For instance, a fundamental difference of this model with the widely used Chilling Hours Model is that the latter does not consider chill accumulation above 7.2 ºC, nor chilling negation, and quantifies the same effect of chilling for any temperature below (or equal to) 7.2 ºC. For areas with warm winters, like Region of Murcia, this can lead to unrealistic characterization regarding chill accumulation. In the absence of significant mitigation measures, global warming is likely to produce an increase in temperatures between 1.5 ◦C and 2 ◦C within the period 2030–2050 compared to pre-industrial annual mean global temperature levels (IPCC, 2021). The rise in temperatures has decreased winter chill accumulation in multiple fruit-growing regions, negatively impacting fruit and nut production (Baldocchi and Wong, 2008; Benmoussa et al., 2020, 2018; Fernandez et al., 2021, 2020a; Luedeling et al., 2009; Rodríguez et al., 2019), even though recent research indicates that lack of chilling can be partially compensated by spring forcing (Kaufmann and Blanke, 2019), which would potentially reduce this negative impact. The lack of accumulated chill negatively affects yield production in multiple ways: delay and lack of sprouting uniformity, decrease in the number of flowers and vegetative buds, flower drop (Agustí, 2010), an unusual shape of fruits (Li et al., 2016; Wert et al., 2007), delay of flowering dates (Li et al., 2016), fewer flowers, longer flowering period and double fruit in cherry (Drogoudi et al., 2020). The Region of Murcia is located in the Mediterranean basin, one of the climate change hotspots (Alpert et al., 2008; Giorgi et al., 2004; Giorgi and Lionello, 2008; Ozturk et al., 2015). Multiple studies have explored possible negative impacts of climate change on fruit and nut production in different areas of the Mediterranean basin provoked by the lack of winter chill such as Tunisia (Benmoussa et al., 2020, 2018), Spain (Egea et al., 2022; Rodríguez et al., 2021) or Greece (Drogoudi et al., 2020). A recent study evaluated the possible negative consequences of future temperature rise on temperate fruit orchards considering the whole Mediterranean basin (Fernandez et al., 2023). Efficient adaptation measures of the orchards will be fundamental to minimising the negative economic impact of global warming. For example, moving crops to higher altitudes within the same region if possible, or providing accurate winter chill accumulation estimates for different locations obtained through spatial interpolation could help farmers in the choice of varieties not only adapted to current climate but also resilient to climate change. The Region of Murcia is one of the main stone fruit producers in Spain: 303.528 t of peaches, 14.382 t of plums and 64.524 t of apricots were produced in 2020 (Centro Regional de Estadística de la Regi´ on de Murcia, 2022). Besides, it contains quite different areas from the agroclimatic point of view: from cold areas in Northern and North-Western locations to warm areas like those close to the coastline. As an illustration of the wide range of agro-climatic requirements of stone fruit trees grown in Murcia, examples of two cultivars of apricot, Japanese plum and peach and their known chilling requirements in chill portions (CP) are given: for apricot "Cebasred" (38 CP), and "Farbaly" (70 CP); for Japanese plum "Black Splendor" (32 CP) and "Songold" (61 CP); for peach "Honey Prima" (25 CP) and "Catherina" (65 CP). These cultivars and all the others present in Murcia are currently grown in their optimum zones in terms of agro-climatic conditions, although some of them could be at risk due to global warming from the point of view of adaptation as they are located in areas where their chilling requirements are currently barely fulfilled. Mapping current and future climate conditions is crucial in agriculture to have an idea of potential shifts in phenological stages of crops (Fraga and Santos, 2021; Grillakis et al., 2022; Parker and Abatzoglou, 2018; Vanalli et al., 2021) to, for instance, make optimal decisions about crop choice (Ad˜ ao et al., 2023) or to design agronomic strategies to mitigate climate change effects (Egea et al., 2022; Freitas et al., 2023). Temperatures or other agroclimatic metrics covering wide areas can come from gridded data (e.g. (Haylock et al., 2008)) or from spatial interpolation using a number of real observations. This latter approach results in lower biases (Liu et al., 2020) and we will refer to it in this study. The number of available weather stations that provide climatic data to characterize agricultural production areas depends on the region. In regions where a dense network of weather stations exists (as is the case of our case study in Region of Murcia, Spain), spatial interpolation techniques are preferred over other options like e.g. gridded weather data to characterize the whole regional territory (Mourtzinis et al., 2017). While there is an abundant number of interpolation techniques, not all are useful to capture the agroclimatic particularities of certain locations, where micro-climates or particular conditions which differ from the expected pattern may happen (Ali et al., 2021; Benavides et al., 2007; Chai et al., 2011; Parmentier et al., 2014). In the context of stone fruit trees and chill accumulation mapping, the choice of the interpolation method can determine the accuracy of the chill accumulation estimates. While there are multiple examples of spatial interpolation of chill accumulation at different scales and regions (Benmoussa et al., 2020; Egea et al., 2022; Fernandez et al., 2023, 2021; Luedeling et al., 2011, 2009; Luedeling and Brown, 2011; Moral et al., 2021; Ribeiro et al., 2024), some of them use the (12-neighbour) kriging method (Luedeling et al., 2011; Luedeling and Brown, 2011) without a prior comparison and have been applied at a global scale. However, recent studies have started exploring other spatial interpolation techniques, such as regression-kriging (Moral et al., 2021) or the spatial interpolation of altitude, an auxiliary variable (Benmoussa et al., 2020; Fernandez et al., 2021). They have been applied to smaller areas, such as Spain, Tunisia or South America, respectively. Taking into consideration that spatial interpolators are not only variable-specific (Yu et al., 2013) but also sensitive to the data structure (Fazeli Sangani et al., 2019; Li and Heap, 2011), it is important not to assume that a method that showed a good predictive power for a region is always suitable for other regions. In order to choose the optimal method, the predictive capability of different techniques should be measured for new data. While for climatic variables (Ali et al., 2021; Benavides et al., 2007; Chai et al., 2011; de Lucas et al., 2023; Mokarram and Mirsoleimani, 2022; Parmentier et al., 2014; Roznik et al., 2019), soil variables (Bishop and McBratney, 2001; Fazeli Sangani et al., 2019) or even other agroclimatic variables (Hodam et al., 2017) comparisons of a wide range of spatial interpolation techniques have been made, a proper comparison between multiple methods is to be made for winter chill accumulation. The inclusion of auxiliary variables in spatial interpolations has been proven useful in climatic studies (Benavides et al., 2007; Benmoussa et al., 2020; Chai et al., 2011; Hengl et al., 2007; Hijmans et al., 2005). Regarding winter chill accumulation spatial interpolation, altitude has been proven useful in the interpolation of minimum and maximum daily temperatures (Parmentier et al., 2014) which is closely related to chill accumulation (Fernandez et al., 2021). In addition, the relationship between altitude and chill accumulation has been proven in field studies (Alburquerque et al., 2008; Cook et al., 2017). In consequence, altitude has been considered in winter chill spatial interpolation, either through the performance of linear regression with altitude as an exploratory variable before the kriging application (Moral et al., 2021) or by directly spatially interpolating altitude and then transforming the altitude map into a chill accumulation map through linear regression (Benmoussa et al., 2020; Fernandez et al., 2020a). Fewer studies have used other auxiliary variables, such as latitude or marine influence, in winter chill spatial interpolations (Moral et al., 2021), even though the direct relationship between latitude and temperature is well-known and the regulatory effect of the sea in both daily and annual extreme temperatures has been documented (Cheval and Constantin, 2019). One of the few studies that consider other explanatory variables apart from altitude, is the one performed by (Moral et al., 2021), in which they evaluated the goodness-of-fit of different linear regressions that explain chill accumulation through altitude, latitude and distance to the sea separately. While altitude explained the highest I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 2 amount of variability, the influence of latitude was also significant, and the inclusion of latitude on the altitude-dependent regression improved the goodness of fit. Comparisons between spatial interpolation methods that use altitude and other exploratory variables and distance-only-based spatial interpolation methods are yet to be made to measure to which extent the consideration of these exploratory variables increases chill accumulation interpolation accuracy. In this study, we review a set of spatial interpolation and prediction techniques to accurately characterize areas at a regional scale from the agroclimatic point of view. We apply the study to the chill accumulation in the Region of Murcia and, given the great importance of stone fruit production in this area. 2. Methods 2.1. Study area, selection of weather stations and gap filling We focus our study in Region de Murcia, located in the South-east of the Iberian Peninsula (~ 38◦0’N, 1◦30 ′ W) with an extension of 11,314 km 2 (Fig. 1). We selected a set of 49 weather stations, 46 of them managed by the Regional Agronomic Information Systems (SIAM) and 3 of them managed by the State Meteorological Agency of Spain (AEMET). The selection criterion was to have at least 15 years of data between 2000 and 2020. Coordinates and elevation of all the 49 selected stations are shown in Table S1. Temperature data were obtained from SIAM and AEMET databases. Hourly mean temperatures were recorded for SIAMs stations. However, daily temperatures were the only metrics available for AEMET stations. In these cases, hourly temperatures were calculated with idealized temperature curves throughout the day (Linvill, 1990) using the chillR package (Luedeling and Fernandez, 2022). Temperature data from AEMET stations of neighbouring provinces were included in order to avoid border effects. In particular, those stations, named using the same coding as AEMET, were 7247X from Alicante province (~ 38◦24 ′ N, 1◦02’W), 7096B from Albacete province (~ 38◦45 N, 1◦40 ′ W) and 6367B from Almería province (~ 37◦24’N 1◦54 ′ W). Almost all the stations had gaps in their temperature reports. For the 49 stations used to calculate chill accumulation, gaps were completed using different methods depending on the gap size. When enough data to accurately estimate daily temperature extremes were available, idealized temperature curves throughout the day were calculated based on latitude and maximum and minimum daily temperatures and they were used to fill hourly gaps. Idealized temperature curves throughout the day were also used to obtain hourly temperatures, needed to calculate chill accumulation. Maximum and minimum daily temperatures are assumed to be known if three hours around the coldest and hottest times of the day are recorded. According to this criterion, maximum and minimum daily temperatures could not be calculated for some days. In these cases, daily temperature extremes were estimated by taking the daily values of the nearest station with available data (as long as the Fig. 1. Climatic stations in the region of Murcia according to their source and their use in this study. I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 3 distance between stations was lower than 25 km, see Figure S1), and then the idealized temperature curves were calculated. The gap-filling process was performed using functions from ChillR package (Luedeling and Fernandez, 2022), as described in (Buerkert et al., 2020). More details about these procedures can be found in (Luedeling, 2018). 2.2. Calculation of chill accumulation Chill accumulation was calculated for each winter season (starting on November 1st and ending on February 28th) between 2000 and 2020, for each of the 49 stations. Chill accumulation was measured in chill portions (CP) using the Dynamic Model. Briefly, the model quantifies chill portions which would be directly related to the amount of an ideal substance formed in the plant allowing to endodormancy release once a certain genotype-dependent value is achieved. Chill accumulates continuously with low temperatures, although the model considers some chilling negation produced by high temperatures. Once a certain amount of chill is accumulated (a chill portion) it cannot further be reduced by negation. The model uses hourly temperatures for its calculation (Fishman et al., 1987a, 1987b). The mean accumulated chill and the Safe Winter Chill (SWC), corresponding to the 10th quantile of a chill accumulation distribution, were computed. The latter can be interpreted as the minimum amount of chilling that growers can expect in 90 % of all years. The authors that introduced the term stated that SWC can be more useful for tree crop growers than mean winter chill, since it incorporates the economic need for an orchard operation to produce good yields in most years (90 %), rather than in an average year (Luedeling et al., 2009). 2.3. Spatial interpolation In this study, a set of spatial interpolation methods with different parameterizations applied to the registered accumulated chill over the studied locations were compared based on their prediction capability. To evaluate predictability, we employed bootstrap cross-validation. This involved partitioning our dataset into a training set and a test set, being the latter set to 10 stations, following the recommendation of using 20 % of data for testing (Marcot and Hanea, 2021). The prediction capability was estimated for the test set, by comparing the chill portions estimated by the model with the observed values through the calculation of the mean squared error. A thousand simulations for each of the spatial interpolation methods were performed through bootstrap, obtaining a mean squared error distribution for each method. In addition, a leave-one-out prediction procedure was performed to identify for which stations winter chill is more difficult to predict. Different spatial interpolation methods (Table 1) were applied to mean and SWC (in portions). A brief description of each of the spatial interpolation methods is presented in the Supplementary Material. Some of the methods in Table 1 can be applied using different parameterizations. For example, IDW contains the power parameter, idp, that relates the influence of the observations with distance. We choose 5 values for this parameter (see Supplementary Material). Likewise, SK and OK contain the so-called “nugget parameter”which is the number of neighbours used to perform the interpolation. Here we explored two values in each case: a high value (n=48) and a low value (n=7), finding no significant differences among them. Besides, all the methods except 3-D Kriging were directly applied to chill portions but also to the residuals of the two regression models used: Portions vs altitude and Portions vs. altitude +latitude (see Table 2 and Fig. 3). All this leads to a total number of instances of 32 considering all the methods, all the tested parameterizations and their application directly to the chill portions data or to the residuals of the selected regression models. To have a compact notation, in this paper we will refer to the methods with the acronyms in Table 1, preceded by "Reg alt" or "Reg alt - lat" if they were applied to the residuals of the altitude or altitude +latitude regression models respectively. In the case of specifying the value of any parameter, it is written in parentheses at the end. For example, "Reg alt +IDW (idp=1.5)" refers to the application of IDW on the residuals of the regression of chill portions versus altitude, with a power parameter equal to 1.5. For temperature data, several studies found that the prediction accuracy increases when spatial interpolation is applied to the residuals of regression rather than directly to the observed temperature data (Bishop and McBratney, 2001; Chai et al., 2011; Hengl et al., 2007; Knotters et al., 1995). By incorporating a multiple linear regression with auxiliary variables, the chill portions’spatial variance explained by other variables can be considered, in addition to the spatial-related variation. Additional exploratory variables related to air temperature and chill accumulation were explored, mainly altitude and distance to the sea, as well as, the interactions between them and with the location. Different multiple linear models with different combinations of the variables: longitude, latitude, altitude, distance to the sea and the interactions between each pair, were generated through the function regsubset from the R package “leaps”and compared based on goodness of fit and simplicity through multiple metrics: adjusted R-squared, Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). Also, distance to the sea squared was included, as in some cases sea influence is concentrated on the first meters of the coast (Cheval and Constantin, 2019). Finally, the spatial interpolation methods previously described were applied. 2.4. Winter chill maps Once the optimum interpolation method was identified, winter chill maps of Murcia were created and compared to the 12-neighbours kriging map, a common spatial interpolation method for winter chill. The resulting winter chill map resolution is limited to the resolution of the explicative variable maps. Elevation maps for Murcia were available on the National Geographic Institute of Spain website. Finally, 25 m resolution mean and SWC maps were produced. 3. Results Mean and safe winter chill were calculated for 49 stations across Murcia. The high diversity of chill accumulation in this region is represented in Fig. 2, which shows some coastal areas with only 54 mean CP like Mazarr´ on (AL62), in contrast with interior and mountainous localities like Caravaca de la Cruz (CR12) where mean chill accumulation reaches more than 75 CP. For both mean and safe winter chill, the inclusion of latitude in the linear regression increases the prediction capability compared to considering only altitude as an independent variable. The inclusion of longitude also improves the statistical metrics AIC and BIC but not as much as with latitude (Table S2). According to the Wilcoxon test, as shown in Table 2, the inclusion of latitude in the linear regression increases prediction accuracy with respect to the model considering only altitude. The prediction capability of each one of the 32 interpolation instances was evaluated through out-of-bag bootstrap. In Fig. 3 we represent the mean squared error distribution of the 1000 bootstrap simulations for each interpolation method for mean winter chill Table 1 Interpolation methods tested in this study. Method Reference Inverse Distance Weighting (IDW) (Shepard, 1968) Radial Basis Functions (RBFs) (Broomhead and Lowe, 1988) Simple kriging (SK) (Matheron, 1967) Ordinary kriging (OK) (Matheron, 1967) Empirical Bayesian Kriging (EBK) (Krivoruchko et al., 2011) 3-D kriging (Goovaerts, 1999) I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 4 (Fig. 3A) and SWC (Fig. 3B), respectively. 3.1. Bootstrap-based results Interpolation methods applied directly are less precise than those applied to the residuals of linear regression including altitude, as shown in Table S3, where the median values of the mean squared error distribution of the 1000 bootstrap simulations for different interpolation instances is presented. These values are presented for both mean and safe winter chill. For each variable, the first column presents the value for the application of the methods directly to chill portions, while the second and third columns present the values corresponding to the application of the methods to the residuals of the altitude and altitude+latitude regression models, respectively. Performing linear regression before spatially interpolating not only improves prediction but decreases performance variability across different sets of testing stations. Mean winter chill is spatially dependent as it is shown in the variogram (Figure S2, top-left). There is also evidence of slight spatial dependency of the altitude dependent regression residuals (Figure S2, middle-left). In consequence, the application of some spatial interpolation methods increases altitude dependent regression prediction accuracy in terms of mean squared error (MSE). In particular, RBF and IDW with a power parameter value (idp) equal to 1 or 1.5 (Table S3) present lower MSE values with statistical differences with respect to regression. On the other hand, altitude+latitude dependent regression residuals do not show any evidence of spatial dependency (Figure S2, bottom-left). None of the spatial interpolation methods applied to the residuals of the altitude+latitude dependent regression statistically outperformed the prediction accuracy of the regression. Table S4 shows the p-values of the statistical tests that compare the mean squared error distribution (calculated through the out-of-bag bootstrap method) of the best spatial interpolation methods for mean winter chill. The methods considered were, for altitude dependent regression residuals, the reference method (i.e., altitude-dependent regression) and those statistically outperforming it: IDW with idp = {1, 1.5} and RBF. For altitude-latitude dependent regression residuals, no method statistically outperformed the altitude +latitude regression, thus this method was the only one considered for the comparison. Table S4 shows that, as stated above, altitude regression-based methods IDW with idp ={1, 1.5} and RBF outperform the altitude dependent regression. Altitude +latitude dependent regression outperforms all the rest of the considered methods but no clear conclusions can be taken on its comparison with altitude regression +RBF (p-value =0.05366). Therefore, both methods will be used for comparison when evaluating mean winter chill. Finally, we explored what resulted in a higher improvement over the altitude regression predictive capability: either adding latitude to the regression or spatially interpolating the residuals with RBF. RBF over the Table 2 Goodness-of-fit and predictive capability of different linear regressions dependent on geographic variables. Linear regression Goodness of fit Predictive capability R 2 R 2 adj AIC BIC Bootstrap distribution median Wilcoxon test p-value Mean winter chill ~altitude 0.714 0.708 256.71 262.32 12.21 2.02⋅10 −09 ~altitude+latitude 0.747 0.735 252.95 260.43 11.40 Safe winter chill ~altitude 0.679 0.672 289.21 294.82 23.95 2.20⋅10 −16 ~altitude+latitude 0.719 0.707 284.78 292.27 21.94 Fig. 2. Mean and safe winter chill accumulation for each climatic station in Murcia, calculated through the Dynamic Model with temperature data from 2000 to 2020. I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 5 residuals of the altitude-dependent regression was slightly more accurate than the altitude+latitude dependent linear regression (Table S3). However, RBF is far more complex and requires much more computational resources than linear regression. As for mean winter chill, altitude+latitude dependent regression provides better goodness of fit than altitude dependent regression for SWC (Table S2). Spatial dependency is also evident in SWC, while altitude-dependent regression residuals are not spatially dependent (Figure S2). Only RBF applied over the altitude dependent regression residuals significantly outperforms altitude dependent regression. None of the spatial interpolation methods statistically increases the accuracy of the altitude+latitude dependent regression. The comparison between RBF applied over altitude regression residuals and altitude+latitude regression concludes that the latter provides a better accuracy (Table S5). In conclusion, the most accurate interpolation method is also one of the simplest: the altitude+latitude dependent linear regression. 3.2. Mapping For mean winter chill, to determine if the slight increase in accuracy offered by the altitude regression+RBF in terms of the mean squared error distribution compensates for the higher requirement of computational resources, two maps of chill accumulation in Murcia were made according to the two different techniques (Fig. 4) for visual comparison. Some differences were found in the warmer areas, near the southern and south-eastern coast but, in general, both maps are very similar. In this study, we will consider Radial Basis Functions of the altitudedependent regression residuals as the optimal method, although altitude+latitude dependent regression is recommended if computational resources are limited as it is much less computationally intensive. The importance of evaluating and selecting the best interpolation method for a given set of data is demonstrated in Fig. 5, where we compare the maps obtained from the first ranked method for each of the Fig. 3. Predictability of each spatial interpolation method used to estimate mean (A) and safe winter chill (B) calculated through out-of-bag bootstrap. I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 6 interpolated metrics (i.e, RBF over the residuals of the altitudedependent regression for mean winter chill and altitude+latitude regression for SWC) with the ones produced by a previously used technique for mapping winter chill: 12-neighbour kriging (Luedeling et al., 2011, 2009; Luedeling and Brown, 2011). While kriging shows a decrease in accumulated winter chill as we move toward the coast and southwards, it does not capture the differences at smaller scales. For example, the noticeable differences between mountainous areas and valleys. The differences in mean squared error, calculated through bootstrap, between both methods, are 19 for mean winter chill and 25 for SWC. As will be explored in the following section, incertitude differences across methods vary importantly between different locations. 3.3. Interpolation accuracy depending on the location: leave-one-out prediction analysis The accuracy of each station’s winter chill accumulation prediction, according to the different spatial interpolation methods, was evaluated through leave-one-out prediction analysis. Fig. 6 represents the distributions of predicted values according to 3 groups of interpolation methods for each of the stations for mean (Fig. 6A) and safe winter chill (Fig. 6B). The real value of winter chill accumulation in each case is represented by a black dot. The groups are formed according to the variables they consider: only distance, altitude and both altitude and latitude. Chill accumulation predictions were less exact for warm locations regardless of the method used, tending to overestimate winter chill (i.e., those with low accumulated chill, Fig. 6). This weakness of the tested interpolation methods is also apparent when we compare real chill accumulation in warm locations with the spatially interpolated maps, as these maps (Fig. 5) do not include values as low as we find in some stations (Fig. 2). The discrepancies between interpolation methods vary highly between stations. Stations in elevated areas are the ones where discrepancies are more significant, especially between interpolation methods that take altitude into consideration and techniques based only on distance. While the latter underestimate chill accumulation, altitudebased methods slightly overestimate it. This discrepancy increases with altitude, as well as the interpolation error. Interpolation methods that include altitude are less precise in elevated stations. Likely, this bias is due to the characteristics of our data set, where localities over 500 m are underrepresented. However, as we have seen through bootstrap analysis, the inclusion of altitude (especially through linear regression) is the main factor to increase winter chill interpolation accuracy. Interpolation methods that consider altitude assign a greater range of values than distance-only methods allowing for more accurate predictions, while distance-only methods tend to underestimate winter chill in cooler locations and vice versa (Fig. 6). Also, there are fewer discrepancies between spatial interpolation methods that include altitude, than between distance-only methods. 4. Discussion The evaluation of spatial interpolation methods’predictive capability through bootstrap showed significant differences between the predictive performance of the methods. While the use and comparison of interpolation methods for different (agro)climatic variables have been already presented within the literature (Ali et al., 2021; Benavides et al., 2007; Chai et al., 2011; Hodam et al., 2017; Mokarram and Mirsoleimani, 2022; Parmentier et al., 2014; Roznik et al., 2019), to the best of our knowledge, it has not been done for winter chill accumulation. According to our analysis, the best-performing method for mean winter chill was RBF applied to the altitude-dependent regression residuals, followed closely by the simpler altitude+latitude regression, which would be recommended if computational resources are scarce. For SWC, the altitude+latitude dependent regression was the best method, with RBF applied to the altitude+latitude regression residuals performing almost identically. These comparisons of methods allowed us to produce Fig. 4. Mean winter chill map of Murcia, calculated through the two best performing methods. Left: altitude-latitude dependent regression. Right: Radial Basis Functions applied over the residuals of the altitude dependent regression. I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 7 a mean winter chill accumulation map with a mean error of 4.3 chill portions (6.4 % of the mean chill accumulation value) and a SWC map with an average error of 6.1 chill portions (7 % of the average SWC value). Considering particular locations, the differences in accuracy between methods are not homogenous across Murcia. Chill accumulation estimates in important locations for agriculture vary widely depending on the method. For example, the prediction error for SWC estimations in Cieza (corresponding to CI42 climatic station), the most important location in the production of peaches in Murcia (Centro Regional de Fig. 5. Winter chill accumulation maps of Murcia, calculated through 12-neighbour kriging (left) and the best ranked method (right) for mean winter chill (A) and SWC (B). I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 8 Estadística de Murcia, 2022), is about 0.8 chill portions according to the altitude-latitude regression (selected as the first-ranked method for SWC after the bootstrap analysis) and 6.6 chill portions when RBF is directly applied to the data. However, as it has been demonstrated through the leave-one-out prediction analysis the differences in the accuracy of the spatial interpolation methods are more notable in higher altitude locations. While most of the agricultural activity takes place between 100 and 400 m of altitude in the Region of Murcia, some cultivars can be located in high-altitude locations, where the choice of the spatial interpolation method is especially important. For example, in Caravaca de la Cruz (corresponding to the CR12 climatic station; 869 m of altitude), the second most important location for plum cultivars in Murcia, the difference in accuracy across methods reaches 11 chill portions for SWC. A similar result was previously observed by (Kilibarda et al., 2014) for the interpolation of daily temperatures, or by (Fernandez et al., 2023) for measuring chill accumulation, obtaining the lowest prediction accuracy at high altitudes. These differences in the predictions among methods can be negligible for some areas, specially taking into account Fig. 6. Chill accumulation (A: Mean winter chill and B: SWC) estimations through different methods in contrast to the real value, represented by a black dot. The Spatial Interpolation Methods are grouped by which information they consider: altitude and latitude (red circles), altitude (blue squares) or only distance (green triangles). I. G´ omez-Ramos et al. European Journal of Agronomy 160 (2024) 127316 9