From intermittent to persistent chill insufficiency in California’s specialty crops

Study area

Our analysis focuses on regions across California where walnuts, pistachios, cherries, and plums are cultivated (Fig. 5). We assessed chill sufficiency separately for each crop, specific to its growing areas. To delineate these areas, we used the 2022 Cropland Data Layer (CDL) from USDA-NASS, which provides land use data at a 30 m resolution⁴². Crop-Specific masks were created by aggregating the CDL data to align with the 4 km resolution of the observational data. Chill sufficiency computation was then generated for each masked region at this 4 km resolution. Additionally, we emphasized crop-growing regions in southern California (south of 37° latitude), where previous studies have reported chill insufficiency risks^16,25. A county-level analysis was also conducted for the six counties with the largest cultivation areas for these crops.

The alternative text for this image may have been generated using AI.

Colored areas in the map represent regions where each crop is grown. The six counties with the largest crop-growing area for each crop are highlighted with bold blue outlines.

Data

Daily maximum (T_max) and minimum (T_min) temperatures from the observational dataset (GRIDMET) were used as the reference for historical observations in this study. GRIDMET provides temperature data at a 4 km spatial resolution, available from 1979 up to the previous day. The dataset is developed by merging temporally rich data from the North American Land Data Assimilation System Phase 2 (NLDAS-2) with the spatially fine resolution (800-m) Parameter-elevation Regressions on Independent Slopes Model (PRISM) data⁴³. It was validated against an extensive network of automated weather stations across the western United States, showing skill comparable to station-based interpolation methods. For consistency with the model forecast period, we selected data from 1981 to 2024.

Model data from the ECMWF, consisting of six-hourly forecasts, were averaged to obtain daily values for the November–December–January–February (NDJF) period. Forecasts were initialized on the first day of each month (November, December, January, and February) from 1981 to 2024 (https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download), including both 1-month lead forecasts (e.g., November initialized in November) and longer lead times up to 4 months depending on the initialization month. Thus, the lead time varied by forecast start date, ranging from near-zero (same-month initialization) to approximately 4 months in advance. Additionally, monthly mean forecasts for the NDJF season were obtained from ECMWF for the same years from one to four-month lead times (https://cds.climate.copernicus.eu/datasets/seasonal-monthly-single-levels?tab=download). We chose ECMWF due to its widely recognized performance, as previous studies have established it as the gold standard for medium-range weather forecasting⁴⁴.

The ECMWF hindcasts are available at a spatial resolution of 1° × 1° (~111 km). However, at this coarse resolution, the model data are limited in their ability to capture temperature variations influenced by topography, proximity to coastlines, and other meso-climatic factors. In contrast, the observational dataset offers finer spatial detail and can represent these small-scale temperature variations more accurately, which is critical for realistic predictions of winter chill. To leverage the large-scale forecast signal while retaining fine-scale variability, we apply the ECMWF anomalies to the 4 km observational climatology using both mean and SD scaling. In this approach, the coarse-scale forecast anomaly is applied across the high-resolution grid, while the local variability from observations is explicitly preserved. This SD-ratio scaling improves chill estimation relative to using raw forecasts or mean-only bias correction. Our aim is not to predict exact daily chill at individual locations, but to capture anomalies in chill sufficiency at monthly-to-seasonal scales for crop-growing regions. This method has been successfully applied in previous studies^16,45,46.

To compute CP, hourly temperatures are required, which were derived from daily maximum (T_max) and minimum (T_min) temperatures. The daily forecasts of T_max and T_min from ECMWF were bias-corrected using the nearest observational grid point and scaled through SD-ratio adjustment to match the observed variance.

In Eqs. (1) and (2), \({T}_{\max D}^{{{\mathrm{BC}}}}\) and \({T}_{\min D}^{{{\mathrm{BC}}}}\) represent the bias-corrected daily maximum and minimum temperatures from ECMWF daily forecasts at a spatial resolution of the observational data, around 4 km for each lead time. The daily climatology \({T}_{\max }^{{{\mathrm{obsDClim}}}}\) and \({T}_{\min }^{{{\mathrm{obsDClim}}}}\) in the observational dataset were derived by averaging daily values for each calendar day over the period 1981–2024. \({T}_{\max }^{{{\mathrm{modDAnom}}}}\) and \({T}_{\min }^{{{\mathrm{modDAnom}}}}\) refer to the daily anomalies of ECMWF’s T_max and T_min, respectively. The daily anomalies were computed by subtracting the corresponding daily climatological means from the ECMWF hindcasts. The terms \({\sigma }_{{T}_{\max D}}^{{{\mathrm{obs}}}}\) and \({\sigma }_{{T}_{\min D}}^{{{\mathrm{obs}}}}\) denote the daily SDs of maximum and minimum temperatures from the observation, while \({\sigma }_{{T}_{\max D}}^{{{\mathrm{mod}}}}\) and \({\sigma }_{{T}_{\min D}}^{{{\mathrm{mod}}}}\) represent the corresponding SDs from the ECMWF forecasts. The daily SD of the corresponding variables was computed for each calendar day over the entire period.

Since the observational data were available at a finer spatial resolution (4 km) compared to the coarser resolution of the ECMWF forecasts (1° × 1°), it was necessary to align their resolutions to perform the arithmetic operations in the equation. To achieve this, we resampled the ECMWF forecasts to match the resolution of observational data by replicating each grid cell across the corresponding observational data grids.

The monthly forecasts of T_max and T_min from ECMWF were then downscaled to daily values using observational data from the nearest grid point to each ECMWF grid cell. Bias and variance scaling (SD-ratio) was applied in this process to preserve the observed interannual variability while retaining the large-scale anomaly signal. This is expressed as:

In Eqs. (3) and (4), \({T}_{\max M}^{{{\mathrm{BC}}}}\) and \({T}_{\min M}^{{{\mathrm{BC}}}}\) refer to the daily maximum and minimum temperatures generated from ECMWF monthly forecasts at a spatial resolution of the observational data, 4 km for each lead time. \({T}_{\max }^{{{\mathrm{modMAnom}}}}\) and \({T}_{\min }^{{{\mathrm{modMAnom}}}}\) are the monthly anomalies of T_max and T_min from ECMWF, calculated by subtracting their respective climatology from the ECMWF hindcasts. The terms \({\sigma }_{{T}_{\max M}}^{{{\mathrm{mod}}}}\) and \({\sigma }_{{T}_{\min M}}^{{{\mathrm{mod}}}}\) denote the monthly SDs of T_max and T_min from the ECMWF forecasts, while \({\sigma }_{{T}_{{{\mathrm{maxM}}}}}^{{{\mathrm{obs}}}}\) and \({\sigma }_{{T}_{\min M}}^{{{\mathrm{obs}}}}\) represent the corresponding SDs from the reference data.

Hourly temperatures needed for chill computation were derived from the daily T_max and T_min, latitude, and sunrise and sunset times. Daytime warming was modeled using a sine curve, while nighttime cooling was modeled with a logarithmic function, following the methods in the ChillR package developed by Luedeling et al⁴⁷. To evaluate whether deriving hourly temperatures from daily T_max and T_min introduces substantial uncertainty in chill accumulation, we computed annual total CP using both methods for three representative stations (Modesto, Davis, and Five Points) using California Irrigation Management System (CIMIS) data (https://cimis.water.ca.gov/). The annual totals from these two approaches were computed using paired t-tests and Wilcoxon tests (implemented in R), which indicated no significant differences, confirming that the derived hourly temperatures provided a reliable estimate of chill accumulation.

Chill quantification

Chill units are a thermal-time metric quantifying exposure to effective chilling temperatures, and dormancy is completed once sufficient chill units have accumulated⁴⁸. It can be measured as one hour at an optimal temperature, as in the Chilling Hours model, or a more complex unit of accumulated chill, as in the Dynamic Model. The total chilling requirement reflects the number of units or hours accumulated within a specified temperature range, which varies by species and cultivars⁴⁹. Chill accumulation can be quantified using the Chilling Hours model, which counts the number of hours with temperatures between 0 °C and 7.2 °C⁵⁰, while ignoring all other temperatures, or by calculating CP using the dynamic model^1,3. The dynamic model is considered particularly well-suited for warming climates, such as that of California, because it represents chill accumulation in a more physiologically based manner than simpler chill-hour models and partially accounts for the reduction of chilling effectiveness at higher temperatures^51,52,53.

We quantified chill accumulation using CP, calculated with the Dynamic_model function from the ChillR package⁴⁷. The dynamic model estimates CP through a two-step process: first, exposure to cold temperatures leads to the formation of an intermediate product, which can either be destroyed by subsequent warm temperatures or strengthened by moderate temperatures. Once a specific quantity of this intermediate product accumulates, it is permanently converted into a CP, unaffected by later temperature changes. The model needs hourly temperature data (in degrees Celsius) over a specified time period and applies experimentally derived constants within an exponential equation to compute the cumulative CP across the entire duration.

We computed CP for the NDJF season total, as well as monthly total for the entire California using observational data, and ECMWF’s monthly and daily hindcasts for the period 1981–2024. For each season, chill accumulation from November and December of a given year was combined with January and February of the following year (e.g., November–December 1981 with January–February 1982). Chill computation of the 2024–2025 season was not possible, as January and February data for 2025 were unavailable at the time of this study.

Analysis of CP

We computed the slope of the interannual total CP in the observational data between 1980 and 2023. For each grid point, the trend slope of total CP was estimated using Sen’s Slope Estimator²⁰. This method computes the slope for each pair of data points as the ratio of the difference in chill values (y) to the difference in time (x). The final Sen’s Slope is the median of all pairwise slopes, providing a robust estimate of trend magnitude that is less sensitive to outliers or non-normal distributions than a simple least squares linear fit. The statistical significance of the trends (1980–2023) was assessed using the autocorrelation-adjusted Mann–Kendall test (Yue–Pilon modification)^21,22. Both the trend slope and modified significance were computed in Python using the pymannkendall package.

A trend was considered statistically significant if the p-value was less than 0.05. We then clipped the crop-growing regions on the gridded chill data, to determine the proportion of each region experiencing significant chill decline. The proportion of each region experiencing a significant decline was calculated by dividing the number of significant grid cells within the region by the total number of grid cells in that region, then multiplying by 100. We also calculated chill for Southern California (south of 37° latitude) and for individual counties by averaging chill values across all grid cells within the latitude band or county boundaries.

We computed the statistical significance of the difference in the coefficient of variation (CV) of CP between the first (1980–1989) and last (2014–2023) decades. The CV was computed separately for each period, and a statistical comparison was conducted using an independent t-test^54,55. The difference was considered significant if the p-value was less than 0.05. Additionally, we assessed the interannual variability by computing the SD of CP using a 10-year moving window for each crop-growing region.

To examine whether CMIP6 projections capture the observed decline in chill accumulation, we calculated CP using T_max and T_min projections from Global Climate Models (GCMs) participating in Coupled Model Intercomparison Project phase 6 (CMIP6) experiments for both historical (1981–2014) and future (2015–2100) periods for the entire California. We selected three GCMS — ACCESS-CM2, EC-Earth3, EC-Earth3-Veg—based on their demonstrated effectiveness for the California region, as determined through process-based and statistical metric⁵⁶. These models were also chosen because high-resolution (~3 km) statistically downscaled projections are available. These CMIP6 GCMs typically operate at coarse spatial resolutions, they do not adequately capture fine-scale geographic variation critical for accurate chill estimation. Therefore, we used downscaled outputs of these selected GCMs generated through the Localized Constructed Analogs, version 2 (LOCA2)⁵⁷, which provide improved spatial resolution (3 km). The downscaled, gridded datasets for the period 2015–2100 are publicly available in NetCDF format through Cal-Adapt (https://cal-adapt.org/data/download/).

For future projections, we used the high-emission scenario, Shared Socioeconomic Pathways (SSP 5–8.5), which assumes continued growth in greenhouse gas emissions throughout the century⁵⁸. Although multiple scenarios are commonly assessed in climate impact studies, we selected SSP 5–8.5 to enable a direct comparison with earlier winter chill assessments that used the highest emissions scenarios available at the time to establish the earliest plausible timeline of chill decline. Consistent with this objective, SSP 5–8.5 is used here to evaluate whether the upper-bound projected decline in winter chill is sufficient to explain the observed historical trends. Extending this analysis to multiple SSPs to assess scenario-dependent uncertainty remains an important direction for future work.

We clipped these CMIP6 projected chill for each crop-growing region in Southern California and computed the average annual chill accumulation for each crop. We compared the number of years with insufficient chill in observed data and historical projections for the same period (1981–2014) to examine whether CMIP6 projections under- or overestimate the observed decline in chill. Also, we examined whether the continuous decline in chill is similar across all counties.

The chill sufficiency thresholds for pistachio, walnut, cherry, and plum were adopted from Pope⁵. For each crop, chilling requirements vary substantially across cultivars, locations, years, and experimental methods, particularly for stone fruits⁵⁹. We use a crop-specific chill threshold defined as the median across all varieties of a given crop reported in Pope⁵. We use the median of CP, since it is less sensitive to extremes and enables consistent regional-scale comparisons in the absence of county-level varietal data. These thresholds are treated as representative benchmarks rather than fixed physiological limits. Consequently, our results are interpreted in terms of increasing frequency and proximity of episodic chill insufficiency driven by interannual variability, rather than the precise timing of persistent threshold crossings inferred from long-term trends. While cultivars with higher-than-median chill requirements may face greater vulnerability and lower-chill cultivars may experience reduced risk, our objective is not to predict cultivar-specific outcomes but to characterize regional-scale chill risk and its temporal evolution under observed conditions, climate projections, and forecasts.

To evaluate interannual predictive performance, we compared how many years in the observational record were correctly captured by ECMWF forecasts. Forecast accuracy was evaluated at two thresholds: predictions were considered correct if the predicted CP were within 10% or within 5% of the observed NDJF totals, since large overestimates could misrepresent chill sufficiency. Although crops may tolerate up to a 20% shortfall, we adopted this conservative threshold to offer more robust forecasts of chill conditions²⁵.

We evaluated the accuracy of chill forecasts derived from bias-corrected ECMWF hindcasts initialized on the first day of each month (November through February), using both daily and monthly mean products. The analysis included forecasts with lead times ranging from 1 to 4 months, depending on initialization month. Forecast skill in predicting chill was assessed at a shorter 1-month lead, with February forecasts first, followed by comparing longer lead time forecasts.

We performed this analysis at the crop-growing region level in Southern California (<37° latitude) by clipping the chill data to each region and calculating regional averages. Also, we performed the analysis at the individual grid level. Additionally, we calculated the ACC for each grid between the ECMWF-derived chill and that computed from observations.

We also assessed the extent to which interannual variability in winter chill accumulation is influenced by the El Niño–Southern Oscillation (ENSO). To do this, we used the ONI from NOAA’s Climate Prediction Center (CPC), which is based on three-month rolling averages of sea surface temperature (SST) anomalies in the Niño 3.4 region (5° N–5° S, 120°–170° W), derived from the Extended Reconstructed SST dataset, version 5 (ERSST.v5) (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). We then computed both lead and lag correlations between the Niño 3.4 index and CP to evaluate the temporal relationship between ENSO phases and chill variability.