The workflow of this study is articulated in four steps (Figure 1). Historical weather series (1984–2018) in the study area were processed to compute agrometeorological indices, which were used to identify the clusters of hazelnut growing seasons sharing similar climatic conditions via multivariate analyses (step 1). A process-based simulation model to estimate the incidence of rotten hazelnuts was developed using available knowledge on Diaporthe spp. (Emmett et al., 1992; Erincik et al., 2003; Anco et al., 2013): the main components of the epidemiological cycle were formalized in sub-models driven by hourly weather variables and biologically meaningful parameters (step 2). The new model was subjected to a global sensitivity analysis to identify the most relevant parameters in modulating the key outputs, also considering their uncertainty under contrasting climatic conditions from step 1 (step 3). The most relevant parameters were then calibrated using the ground truth data of rotten incidence from post-harvest analysis and phenological observations collected in 22 orchards; the model evaluation was carried out on independent and additional 321 datasets, before running spatially distributed simulations over the whole hazelnut producing area in Turkey (step 4).

The input weather data for the modeling activities obtained from the National Aeronautics and Space Administration (NASA) Langley Research Center (LaRC) Prediction of Worldwide Energy Resource (POWER) Project funded through the NASA Earth Science/Applied Science Program, which provides daily meteorological variables at 0.5° × 0.5° resolution grid. We used here maximum and minimum air temperature (°C), dew point temperature (°C), relative humidity (%), and average wind speed (m s^–1). Hourly air temperature was estimated from daily maximum and minimum air temperature, according to Campbell (1985), whereas air relative humidity was derived according to Bregaglio et al. (2010), based on the models proposed by Linacre (1992), Allen and FAO (1998), and Hahn et al. (1998). Hourly leaf wetness was estimated from hourly air temperature, dew point temperature, wind and relative humidity, according to Kim et al. (2002). In total, 27 NASA-POWER grid cells have been selected to cover the main hazelnut producing regions of Turkey (Figure 2).

The hazelnut rotten simulation model was calibrated using phenological observations and post-harvest rotten incidence data collected in 22 hazelnut orchards located in the municipalities of Samsun, Ordu, Giresun, and Trabzon (Eastern Black Sea area), and Duzce, Sakarya, and Zonguldak (Western Black Sea area) (Figure 2). Phenological observations were collected weekly in these orchards in the 2018 and 2019 growing seasons, whereas the post-harvest analyses of rotten incidence were carried out in the period of 2016–2019 (4 years). The model evaluation was carried out on 321 additional locations where post-harvest analyses on hazelnuts samples were available in the same reference period (2016–2019, as shown in Arciuolo et al., 2020, this special issue). In all datasets, hazelnuts were collected at full ripening (phase R13, as shown in Supplementary Table 2) on 100 trees per orchard (30 hazelnuts per tree), when average kernel humidity was about 10%, to obtain approximately 15 kg of hazelnuts sample^–1 (about 3,000 hazelnuts). Hazelnut fruits were dried on pallets to facilitate the husk removal and mechanical dehusking; after a second drying period to ease the cracking procedure, they were manually shelled and observed for defects after cutting kernels in two halves. The percentage incidence of rotten hazelnuts with visible and invisible defects was determined in laboratory and used as reference data to evaluate model performances.

The 22 orchards where calibration activities were performed are located in 13 out of the 27 NASA grid cells covering the whole hazelnut producing area in Turkey. Therefore weather data from these 13 grid cells were used for model calibration. The model was then evaluated comparing simulation results obtained on all 27 grid cells with median rotten hazelnuts incidence from the additional 321 sampling locations, after averaging them based on the NASA-POWER grid cell they fall in.

Differently from other plant diseases whose etiological agent is well known, several genera (Diaporthe, Alternaria, Cladosporium, Fusarium, and Colletotrichum) are associated with rotten hazelnut defects. Diaporthe spp. emerged as candidate pathogens in the Caucasian region (Battilani et al., 2018). Based on this work, the rotten simulation model presented here is composed by two modules: a process-based model to simulate the epidemiological processes of the Diaporthe spp. cycle, coupled with the reproduction of the susceptibility of hazelnuts to fungal infections as modulated by their phenological development. The large uncertainty associated with the causal agent and the heterogeneity of Diaporthe strains isolated on rotten hazelnuts in Turkey (Arciuolo et al., 2020) led us to develop a hybrid approach, where the indicators of the suitability of weather conditions to generic fungal pathogens are derived from long-term simulations (climatic norm, 30 years over whole Turkey, 1988–2018).

The workflow of the rotten hazelnut simulation model is presented in Figure 3. Model simulations start on October 1, defined as the starting date of the hazelnut growing season. The simulation of the hazelnut reproductive phases is performed according to Bregaglio et al. (2016, 2020, 2021), this special issue. Eight phenological phases were simulated for female reproductive development, from flowering to nut dropping. The simulation of the suitability of weather conditions to rotten increase starts with the formation of pycnidia, i.e., suitable structures for the overwintering of the pathogen (Arciuolo et al., 2021), as driven by air temperature and relative humidity. When pycnidia are formed and mature, they produce cirrhi where conidia develop, i.e., the asexual form of the pathogen (Pscheidt and Pearson, 1991). The model considers that conidia can spread based on the rainfall intensity and duration. Once the plant susceptibility period is reached, the flowers are receptive for the fungal spores. Fungal infection is simulated when conidia are spread, with suitable leaf wetness duration and temperature. The simulated rotten incidence (%) is defined by the following variables: the number of hours suitable for infection events in the current season, the cumulated rainfall in January, and the number of hours suitable for conidia spread in the previous season, to consider the carry over effect of the inoculum load from 1 year to the next.

At each time step, when hourly air relative humidity is higher than a threshold (RHthresholdPC) and temperature is within the minimum and maximum for pycnidia formation (TmaxPF, TminPF), hourly temperature is cumulated. Once the cumulated hourly temperature (CumT, °C) reaches a given threshold (ThresholdPF), pycnidia start to form (PycnidiaFormation, 0-1, Eq. 1).

Conidia are spread (ConidiaSpread, 0-1, Eq. 2) when both cumulated rainfall (CumHourlyRain, mm) and consecutive rainy hours (RainyHours, h) exceed a threshold (RainStartCS; HoursStartCS).

The number of suitable hours for conidia spread (CumCS) is used as a proxy of the inoculum load for the following year and to modulate the increase of rotten incidence (%, IncreaseR, Eq. 3) associated with each infection event. IncreaseR depends on the maximum rotten incidence increase (IP, %), modulated by the normalization of CumCS according to the median of simulated hours when conidia were spread in the period 1984–2018 (CS¯, Eq. 3).

The plant susceptibility (PS, 0-1, Eq. 4) is activated when female flowers are receptive for fungal spores (DVSstart). The codes and description of the phenological phases are reported in Supplementary Table 2. The optimum susceptibility and the end of the susceptible period are driven by the parameters DVSopt and DVSend.

The infection events are simulated as a function of leaf wetness, which is cumulated (CumLW) when hourly air temperature (T, °C) is in the range between the minimum (TminCI, °C) and the maximum (TmaxCI, °C) for infection. Once CumLW and temperature are conducive, a moisture function [f(M), 0-1, Eq. 5] triggers the infection events (Infections, 0-1, Eq. 7). A dry period, i.e., no leaf wetness, terminates the infection when its duration exceeds a threshold (LW50CI):

Where ToptCI is the optimum temperature for infection. Each infection event contributes increasing rotten incidence (Eq. 8). The parameter RottenCoefficient (Eq. 9) is modulated by the inoculum load from the previous year (IncreaseR, Eq. 3) and by the cumulated rainfall in January (CumR), normalized considering the median of the period 1984–2018 (CumR¯):

The rotten model was developed as a BioMA component in Microsoft C#, following the guidelines of the Diseases components (Bregaglio and Donatelli, 2015; Valeriano et al., 2021; Wang et al., 2021).

The seven most relevant parameters explaining the variability of the outputs from the sensitivity analysis assessment were adjusted via automatic calibration, moving their values within their biological ranges (Table 2). The remaining 11 parameters were set as their default values (Table 2). Ground-truth data from the 22 orchards where weekly phenological observations in 2018–2019 and rotten incidence in 2016–2019 were available and were used for model calibration. The rotten incidence model was coupled with a multi-start downhill simplex algorithm (Nelder and Mead, 1965) to perform automatic calibration, setting a weighted root mean square error (RMSE, 0.5) and Pearson’s correlation coefficient (r, 0.5) as objective function to maximize both accuracy and correlations with reference data. We used 10 simplexes and 1,000 iterations, setting 0.001 as the tolerance value (Giuliani et al., 2019). After calibration, the model was applied on an independent dataset of 321 orchards, whose data were aggregated according to the respective NASA-POWER grid cell used as the source of weather data. Model performances were then assessed at municipality level, by taking the median of model results from all the NASA-POWER grid cells of interest. Model accuracy was assessed via mean absolute error (MAE, %), RMSE (%), Pearson’s r, and the coefficient of determination (R²).

Model performances in calibration are reported in Table 4 as aggregated to the municipality level, considering the location of the 22 orchards where weekly phenological observations and rotten incidence data were available (Figure 2). In the Western Black Sea region, RMSE in predicting rotten incidence was always below 2.00%, and MAE ranged between 0.50% in Duzce and 1.40% in Zonguldak. Pearson’s correlation coefficient was higher than 0.50 in all municipalities, even if significant values were reached only in Sakarya, due to the low number of orchards in each municipality. In the Eastern Black Sea region, model accuracy was slightly lower, with RMSE in predicting rotten incidence between 1.26% in Trabzon and 2.55% in Samsun, the latter municipality leading to a MAE higher than 2.00%. Trabzon was the only municipality where Pearson’s correlation coefficient was significant. However, when considering all orchards in the same region, Pearson’s r was significant both in Western (r = 0.72) and Eastern (r = 0.45) region, as well as when evaluation metrics were computed on the whole calibration dataset. Model performances in calibration at the country level denoted a good model accuracy (RMSE = 1.26% and MAE = 0.99%) and a significant correlation (r = 0.58) with field observations.

The model accuracy in reproducing reference data coming from an independent dataset of 321 additional orchards are presented in Figure 7. The histogram analysis highlights that the model correctly reproduced a very high rotten incidence in 2016 in the Western Black Sea region, even if the exceptional rotten incidence in Duzce (>13.00%) was underestimated. Coherently with field observations, the model simulated lower rotten incidence in the period of 2017–2019 in all Western municipalities, with average model errors lower than 2.00% in all cases. In the Eastern Black Sea region, the model was able to reproduce an overall higher rotten incidence than in the Western region, as well as to match the higher observed rotten incidence in 2016–2018, especially in Giresun, Ordu, and Trabzon, whereas in Samsun the model underestimated the reference data. Simulations in 2019 highlighted low rotten incidence in all the municipalities of the Eastern region, in agreement with field data. The evaluation metrics computed on the validation dataset denoted slightly lower performances than in calibration, with RMSE comprised between 2.04% in the Eastern and 2.94% in the Western Black Sea region, and MAE lower than 2.00% in all cases. Pearson’s correlation coefficient computed at municipality level ranged between 0.46 in the Eastern and 0.63 in the Western Black Sea, with corresponding R² of 0.26 at country level.