Impacts of large-scale teleconnection indices on chill accumulation for specialty crops in California

• Large-scale teleconnection showed significant correlations with chill accumulations in California. • The individual teleconnection indices explained up to 39% variability of chill


Introduction
California is the largest and most diverse agricultural state in the United States (US), with an agricultural market value of over $50 billion (USDA NASS, 2019; CDFA, California agricultural statistics review [2019][2020]. The state is the leading, and in some cases, the sole producer of specialty crops in the US, that include fruits, vegetables, tree nuts, dried fruits, horticulture and nursery crops, including floriculture (USDA, Agricultural Marketing Service, 2004). More than 30% of the US cultivation area of specialty crops is in California (USDA NASS, 2019), where more than 65% of US-grown fruit and nuts, including more than 90% of US-grown almonds, avocados, grapes, lemons, mandarins, nectarines, olives, pistachios, plums, and walnuts are produced (CDFA, California agricultural statistics review 2018-2019).
Despite being highly productive, California's agriculture sector is impacted by climate change in numerous ways, from lengthening the frost-free season and increasing heat stress and pest pressure, to shortening the time to crop maturity and reducing soil moisture and irrigation water availability (Pathak et al., 2018). Among the changes in climate, rising temperature is a significant one. The minimum temperatures are increasing at a higher rate than the maximum temperatures for many locations across the state (Ladochy et al., 2011), and especially for the Sacramento-Delta region and South Interior region (Cordero et al., 2011). For California's agriculture sector, one consequence of warming minimum temperatures is lower chill accumulation in winter. Although warming winters may be a potential boon for cool-climate locations (e.g., Gornall et al., 2010), it may be detrimental to crop production in warm-temperate and subtropical locations, like California, where perennial fruit and nut crops rely on cool winter temperatures and enough chill accumulation for their development (Luedeling et al., 2009). Sufficient chill accumulation is essential for the optimal growth and development of many perennial specialty crops. During cool months, woody perennial crops enter endodormancy, a phenological phase during which the plant is hardened to cold temperatures. In response to unfavorable external environmental conditions during this period, internal hormonal regulation prevents bud development, and allows perennials to "track" their exposure to cool temperatures. After a certain, cropspecific amount of time exposed to cool temperatures, also known as chill accumulation, the plant's internal accounting provides a seasonal cue that allows the plant to exit dormancy and continue its development through bloom and fruiting (Saure, 1985). In the absence of sufficiently cool winter temperatures that allow for adequate chill accumulation, crops may experience delayed or uneven budburst, reduced flower quality and inhibited pollination, stunted fruit development, and reduced fruit quality and yield (Atkinson et al., 2013;Luedeling et al., 2009;Saure, 1985).
The interannual variability of winter temperatures, and therefore in chill accumulation, is driven in part by the changes in oceanic and atmospheric circulations, known as teleconnection indices. Teleconnections have been linked to variations in seasonal rainfall, temperature and jet stream anomalies in much of the U.S. (NOAA/National Weather Service, 2008). Their effects are manifested as regional droughts (Rajagopalan et al., 2000), wildfires (Fasullo et al., 2018), floods (Villarini et al., 2013), and crop yield variations (Pathak et al., 2012). In California, precipitation is found to be modulated by teleconnection with significantly increased rainfall during strong El Niño events (Hoell et al., 2016;Bayr et al., 2019). Benson et al. (2003) revealed the Pacific Decadal Oscillation (PDO) maxima had a negative impact on the hydrologic balance of the Sierra Nevada in California and Nevada. Liu et al. (2017) found the Pacific North American pattern (PNA) is strongly correlated with fluctuating drought conditions in the western United States (including California) over the past millennium.
Previous studies have mainly focused on the impacts of teleconnections on individual climate variables, such as precipitation and temperature. However, little attention has been given to the impact of climate variability on integrated climate indices relevant to agriculture, such as chill accumulation. Integrated climate indices are non-linear derivations of the basic climate variables. They can translate raw data into meaningful agricultural indicators, which growers can utilize for immediate and long-term strategic decisions. A recent study also indicates the projections of integrated climate metrics are more useful to farmers than projections of basic climatic metrics for agriculture adaptation (Jagannathan, 2019). Therefore, the study of the impact of climate variability (in the form of teleconnections) on decision-relevant integrated climate metrics fills both a research gap and provides a more useful output for decision making in agriculture.
The objectives of this study are to first assess observed trends in winter chill accumulation in California for two chill accumulation models, and then to evaluate the influence of large-scale teleconnection indices on winter chill accumulation for three different teleconnection indices and their respective pattern phases. The results can be used to develop decision support tools to manage risks for agricultural producers, crop consultants, and policymakers in California, and the methodology can be adapted for other regions to similarly support agricultural decisionmaking. The paper proceeds as follows: Section 2 presents our study area, the temperature dataset and teleconnection indices that we used. Section 3 provides the details on chill accumulation models, method to classify teleconnection modes, and statistical approaches. Section 4 illustrates the results, observed trends in chill accumulation over the last 41 years , and impacts of three teleconnection indices on chill accumulation in California. Section 5 discusses the results, feasible management and adaptation strategies, and the broader impacts on sustainable agriculture, while Section 6 offers concluding remarks.

Study area
California's unique Mediterranean climate, large-scale water transference, and associated irrigation systems make the state one of the most productive fruit and nut growing regions in the US. The state accounts for more than half of specialty crop production nationwide (CDFA, California agricultural statistics review 2019-2020) and over 13% of the nation's total agricultural value (CFDA, 2020). For many specialty crops, such as almonds, wine grapes, strawberries, pistachios, leaf lettuce, and walnuts, California's share approaches or exceeds 90% of US production (Starrs and Goin, 2010; CDFA, California agricultural statistics review 2019-2020). The cultivation of almonds, pistachios, and walnuts are largely concentrated within the Sacramento Valley and San Joaquin Valley, which collectively comprise California's Central Valley. The cultivation of grapes falls within both Central Valley and the arable Central Coast. Strawberries are mainly cultivated in Sacramento Valley and Bay Area, while lettuces are mostly planted in San Joaquin Valley as well as Imperial county (USDA NASS, 2021). Fig. 1 illustrates the basic land cover types based on the National Land Cover Database (NLCD) land cover 2016 product (Yang et al., 2018). According to the 2017 Census of Agriculture (USDA NASS, 2019), about 25% of the state's approximately 100 million acre land area is used for agriculture, of which 47% is permanent pasture and rangeland, 39% is cropland, and about 8% is woodland. Considering the largest specialty crop production centers are the Sacramento Valley, San Joaquin Valley and Central Coast (highlighted in Fig. 1), this study focuses on the chill accumulation analysis within these three regions.

Temperature data
In this study, the gridded temperature data from gridMET were used. gridMET was generated by blending spatially-rich data from the N. Zhang, T.B. Pathak, L.E. Parker et al. Science of the Total Environment 791 (2021) 148025 Parameter-elevation Regressions on Independent Slopes Model (PRISM) with temporally-rich data from the North American Land Data Assimilation System Phase 2 (NLDAS-2) using climatically-aided interpolation (Abatzoglou, 2013). It provides daily updated surface meteorological data covering the contiguous US at 4-km spatial resolution from 1979 to present, which include maximum and minimum temperature, precipitation accumulation, downward surface shortwave radiation, wind velocity, and humidity. Considering the chilling season in California is from November to February, the maximum and minimum daily temperature from Nov.1, 1979 to Feb. 29, 2020 were extracted in this study for calculating chill accumulation.

Teleconnection indices
Previous studies have revealed that the climate variability in the Western US has been influenced by the teleconnections spanning decadal (PDO), interannual (El Niño-Southern Oscillation, ENSO), and intramonthly (PNA) timescales (Abatzoglou, 2011;Bayr et al., 2019;Liu et al., 2017). In this study, three teleconnection indices, including the Oceanic Nino Index (ONI), PDO and PNA were investigated for their impact on the chill accumulation in California. All three indices were extracted from 1979 to 2019, which is consistent with the period of the chill accumulation.
ONI is NOAA's primary indicator for monitoring ENSO (El Niño and La Niña) events. The ONI tracks the sea surface temperatures (SST) in the east-central tropical Pacific between 120°-170°W and 5°N-5°S (Rasmusson and Carpenter, 1982), which is also known as the Niño 3.4 region. To calculate the ONI at certain month, the SST from that month, its previous and following months were required, based on which the 3-month running mean was calculated and compared to a 30-year climatological average. The anomalies or the departure from the climatological mean was the ONI value for the centered month (Trenberth, 1997). The monthly ONI index was obtained from National Weather Service, Climate Prediction Center, NOAA (https://origin.cpc. ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). The operational definition used by NOAA to classify El Niño or La Niña events is the ONI anomalies that exceed ±0.5°C for at least five consecutive months (Trenberth and Stepaniak, 2001). In this study, we adopted a simplified version to classify the positive, negative and neutral modes using ONI, the details of which were provided in the Methods section.
PDO is an El Niño-like oscillatory pattern of Pacific climate variability. It is defined as the leading empirical orthogonal function (EOF) of monthly SST anomalies over the North Pacific basin (north of 20°N) (Mantua, 1999). PDO has considerable influence on the climate in Pacific Basin and over North America, and its signal is most energetic during the winter and spring months (Mantua, 2001). When SSTs are anomalously warm in the interior North Pacific and cool along the Pacific Coast, or sea level pressures are above average over the North Pacific, PDO has a negative value. When the Pacific anomaly patterns are reversed, PDO has a positive value (Mantua, 1999;NOAA/NCDC, 2021). Although the spatial pattern of PDO resembles that of ENSO (Zhang et al., 1997), PDO requires a long data record to compute (e.g., more than 50 years), and it presents the climate variation from interannual-to-interdecadal time scale (Mantua and Hare, 2002). In this study, the PDO dataset generated by NOAA's National Centers for Environmental Information (NCEI) using NOAA's extended reconstruction of SSTs was used (https://www.ncdc.noaa.gov/teleconnections/ pdo/). The NCEI PDO index closely follows the Mantua PDO index.
PNA is the second leading mode of a Rotated Empirical Orthogonal Function (REOF) analysis of monthly mean 500-hPa height anomalies over 0-90°N latitude (Barnston and Livezey, 1987). PNA intrinsically operates on intramonthly timescales, with a life cycle of about two weeks (Feldstein, 2002). The positive phase of PNA is associated with above-average temperatures over western Canada and the extreme western United States, and below-average temperatures across the south-central and southeastern US (NOAA/National Weather Service, 2012). PNA has the largest variability during the cold season, while its signal is least apparent during summer (van den Dool et al., 2000). In this study, the monthly PNA index was obtained from the National Weather Service, Climate Prediction Center, NOAA (https://www.ncdc. noaa.gov/teleconnections/pna/).

Methods
Fig. 2 presents the flowchart of this study. A detailed description of the approach at each step is introduced in this section.

Calculation of chill accumulation
In California, chill accumulation measurement is traditionally from November to February (Luedeling et al., 2009). Multiple models are available for estimating winter chill accumulation, including the Chill Hours Model, Utah Model, Dynamic Model, and Positive Utah Model. Considering the Chill Hours Model is one of the oldest, simplest and most commonly used models (Weinberger, 1950;Jagannathan et al., 2020), and the Dynamic Model is more accurate than other models and widely adopted in California (Luedeling, 2012), these two models were used in this study for the calculation of chill accumulation.

Chill Hours Model
The Chill Hours Model was initially developed in the 1930s. The simplest form of this model accumulates the number of hours with temperatures <7.2°C (Weinberger, 1950) and weighs all temperatures below this threshold equally when estimating the physiological influence on the plant. A modified version only accumulates temperatures between 0°C and 7.2°C (T 7.2 ), assuming temperatures below 0°C or above 7.2°C to have little or no influence on the plant's physiology with respect to chill accumulation (Luedeling et al., 2009). In this study, the modified Chill Hours Model was used to calculate chill hours (CH).

Dynamic Model
The Dynamic Model was originally developed for the warm winters in Israel (Fishman et al., 1987a,b). It allows for warm temperatures to cancel the effects of cool temperatures and accumulates chill in a twostep process that prevents negative chill totals. The first step in this process develops an intermediary chill product that can be built through exposure to cool temperatures and can be destroyed through exposure to warm temperatures. This product transforms into an irreversible "chill portion" (CP) when it reaches a chill threshold. In subtropical and mild temperate climates such as those in California, the Dynamic Model has been to be more accurate than other chill models (Luedeling, 2012).
The R package ('chillR') developed by Luedeling et al. (2013) was used to run the Chill Hours and Dynamic Models in this study. For each of these models, accumulated chill was calculated using hourly temperature downscaled from the daily gridMET temperature time series. The idealized temperature curves proposed by Linvill (1990) were used to convert the daily temperature to hourly data. The Linvill equations use a sine curve to fit daytime temperatures and a logarithmic decline curve for nighttime cooling. The required inputs include daily maximum and minimum temperatures, the latitude of the location to calculate the sunset and sunrise hours, as well as daylength, which we calculated based on the method proposed by Forsythe (1995). The chill units were accumulated from November 1 to January 31 of the next year. In this study, the year of the chill accumulation was defined as the starting year of the chill accumulation period. For example, the chill accumulation from Nov. 2019 to Jan. 2020 is the chill accumulation for year 2019.

Sensitivity analysis of the Chill Model
In this study, the sensitivity analysis was performed on the Chill Hour Model to investigate the sensitivity of chill accumulation in response to the changes in the upper and lower temperature thresholds used in the Chill Model. By comparing different chill accumulation models, Luedeling et al. (2009) found that the Chill Hours Model was most sensitive to rising temperatures compared with the Utah Model and the Dynamic Model. Therefore, the Chill Hour Model was used here as an example to draw the boundary of maximum uncertainty of chill variations. Two experiments were conducted: in first case, a fixed upper threshold (7.2°C) was used, and the lower threshold varied from −3°C to 3°C with 0.5°C intervals. In the second case, a fixed lower threshold (0°C) was adopted, and the upper thresholds varied from 4.2°C to 10.2°C with 0.5°C intervals.

Classify different teleconnection modes
Based on correlation analysis, a simplified rule was adopted to classify the 41 years (1979-2019) into the positive and negative modes using the month that presents the highest correlation between chill accumulation and each teleconnection index. For ONI, a neutral mode is also assessed. For example, when using ONI, the positive mode is the years with ONI > 0.5 at specific month; the negative mode is the years with ONI < −0.5 at a specific month, and neutral mode is the years with ONI values between −0.5°C and +0.5°C (−0.5°≤ ONI ≤ 0.5°C) at a specific month. For both PDO and PNA, the value of zero is used as the cutoff of positive and negative modes.

Statistical analysis
We first calculated the climatological mean of chill accumulation, either for chill hours (CH mean ) or chill portions (CP mean ) by averaging the chill accumulation over 41 years. We then used the non-parametric Mann-Kendall (MK) statistical test (Mann, 1945;Kendall, 1975) to analyze the monotonic trends in the chill accumulation over the 41 years. The null hypothesis of MK test here is that there are no trends in the chill accumulation series, while the alternative hypothesis there is either consistently increasing or decreasing trends in the chill accumulation series at significant level of 0.05. If the chill accumulation series presents significant trend, the slope coefficient in the linear regression model is used to describe the magnitude of the significant trend.
The Pearson correlation coefficient (r) is used in this study to examine the correlation between chill accumulation and teleconnection indices (e.g., ONI, PDO and PNA) at each month. This serves to identify the most influential month of teleconnection indices on chill variation. ONI calculation is based on a 3-month moving window meaning that at the beginning of the chill accumulation period (November 1), the ONI value is only updated to September. Therefore, the monthly correlation between chill accumulation and ONI can only be calculated from January to September. To be consistent with ONI, the correlation between chill accumulation and the other two teleconnection indices (PNA and PDO) are also calculated from January to September. Finally, the paired-sample t-test is adopted to test whether there is a statistically significant difference between two datasets (e.g., averaged chill accumulation at different modes in this study). The null hypothesis of the paired-sample t-test is that the difference between two datasets, (e.g., xy), comes from a normal distribution with a mean equal to zero and unknown variance. show the spatial distribution of CH mean and CP mean in California. It is common for both CH mean and CP mean that the chill accumulation in California is decreasing from north to south. The CH mean of focused regions (Central Coast, Sacramento Valley and San Joaquin Valley) ranges from 147.9 to 1385.3, while CP mean of the focused region ranges from 25.6 to 66.1. Fig. 3(c) and (d) present the magnitude of the linear trend of CH and CP over the last 41 years. It is true for both CH and CP that most of California presents a decreasing trend of chill accumulation, especially for the focused region. After the MK test, a statistically significant decreasing trend of CH is observed over most of the Central Valley (Fig. 3e), while for CP only the northern San Joaquin Valley presents a significantly decreasing trend (Fig. 3f). Over the focused region, the average significant trend of CH is about −4.98/year, and the average significant CP trend is −0.15/year over the 1979-2019 period.

Climatological mean and trend of chill accumulation
It is important to note that the climatological mean and trends in chill accumulation were based on the default upper and lower thresholds used in the literature. Sensitivity analysis revealed that chill accumulation is highly sensitive to the changes in the lower thresholds of Chill Hours Model. Using a fixed lower threshold of 0°C, most of California is experiencing a decreasing trend of chill accumulation (Fig. 3). The average chill hours using a 0°C lower threshold is about 818 CH, while it drops to 598 CH when using a threshold of 3°C. In other words, the 3°C increase of the lower threshold for chill accumulation corresponds to a 27% drop of CH. The increase of the upper threshold may increase the average chill hours by up to 1300. However, the Chill Hours Model does not account for the reversibility of the chill process when temperatures are high (Luedeling et al., 2009). This also agrees with Erez et al. (1979) that the Chill Hours Model did not consider the chill negation effect, in which case the higher temperatures produce negative effects and result in negation of a certain amount of chill units. However, the chill negation is considered in the Dynamic Model, making it more appropriate for chill accumulation study under climate change as noted by Luedeling et al. (2009). More details on the sensitivity analysis,  Table 1).

Correlation between chill accumulation and teleconnection indices
The Pearson correlation (r) between chill accumulation (CH and CP) and three teleconnection indices, including ONI, PDO and PNA, were calculated from January to September. Fig. 4 presents the spatial pattern of the maximum r among the 9 months for CH and CP. When ONI was used, less than 10% of the focused region (Central Coast, Sacramento Valley and San Joaquin Valley) present a significant correlation with CH (Fig. 4a), while about 21% of the focused region show a significant positive correlation with CP, especially at the west and south of San Joaquin Valley (Fig. 4d). ONI explains between 10% and 25% of CP accumulation variability over focused region, with an average of 11%.
PDO presented the strongest correlation with chill accumulation among the three teleconnection indices. PDO explains on average 15% of CH variability in the focused region, ranging between 10% and 30%. The positive correlation between PDO and CH is observed in the Coastal Ranges and Sierra Nevada foothills surrounding the Central Valley ( Fig. 4b), taking about 18% area of the focused region. On the other hand, PDO explains on average of 16% of CP variability, ranging between 10% and 39% over the focused region. Most of the Central Coast and Sacramento Valley (about 43% area of the focused region) are dominated by the significant positive correlation between PDO and CP (Fig. 4e).
PNA shows a roughly opposite pattern to PDO with the most widely spread negative correlation with CH and CP in the focused region ( Fig. 4c and f). Significant negative correlation between PNA and CH (CP) covers approximately 46% (85%) of the focused region. Over focused region, PNA explains between 10% and 29% variability of CH, with an average of 12%; While it explains between 10% and 26% variability of CP, with an average of 15%. Fig. 5 further revealed the corresponding month when the maximum r is achieved. The dominant month when the maximum r is achieved using either CH or CP is September for ONI (Fig. 5a, d), January for PDO (Fig. 5b, e) and September for PNA ( Fig. 5c and f).

Averaged chill accumulation at different teleconnection modes
The teleconnection indices at the months when the maximum r is achieved (Fig. 5) are used to group the 41-years chill accumulation into teleconnection modes (positive, negative and/or neutral) based on the approach introduced in Section 3.2. Using ONI in September, 10 years are identified as the ONI positive mode, 8 years are grouped to the negative mode and the remaining 23 years are classified as neutral. Using PDO in January, 23 years are classified as positive mode and 18 years are grouped to negative mode. September PNA yields 21 years classified as the positive mode and 20 years as the negative mode.
Figs. 6, 7 and 8 present the spatial patterns of the difference between averaged chill at each mode and the climatological mean of chill. The spatial pattern of the averaged chill at each mode is shown in Supplementary Figs. 2, 4 and 6, while the significant differences that passed the paired t-test are shown in Supplementary Figs. 3, 5 and 7. Fig. 6 presents the differences between averaged chill (CH/CP) at each mode of ONI and the chill climatological mean (CH mean /CP mean ). The spatial pattern of the averaged chill at each mode of ONI is shown in the Supplementary Fig. 2, while the statistically significant differences in Fig. 6 that passed the paired t-test are shown in Supplementary Fig. 3. In general, there is no apparent difference of mode-averaged CH/CP among the three modes of ONI ( Supplementary Fig. 2). The modeaveraged CH is generally higher than CH mean in the focused region during both ONI positive (Fig. 6a) and negative (Fig. 6b) mode, while during ONI neutral mode, the mode-averaged CH seems to be slightly lower in the focused region than CH mean (Fig. 6c). However, most of these differences are not statistically significant ( Supplementary Fig. 3a-c), and only 16% area of the focused region show above-average CH accumulation during the ONI positive mode, which mainly lies in the south of Central Coast and San Joaquin Valley ( Supplementary Fig. 3a). In contrast, the mode-averaged CP is slightly lower than the climatological mean in the focused region during both ONI negative (Fig. 6e) and Fig. 4. The maximum correlation (r) from January to September calculated between chill accumulation (CH and CP) and three teleconnection indices (ONI, PDO and PNA). The correlation showed here is statistically significant at p < 0.05. neutral ( Fig. 7f) modes, although the differences are not significant ( Supplementary Fig. 3e-f). However, at ONI positive mode, the focused region is dominated by the positive difference (Fig. 6d), and the significant above-average CP is observed about 57% of the focused region, and mainly in the San Joaquin Valley (Supplementary Fig. 2d). The average significant difference over the focused region is about +3.6 CP, which indicates the focused region may expect 3.6 higher CP during ONI positive mode than the CP mean . Fig. 7 present the difference between averaged chill (CH/CP) at each mode of PNA and the chill climatological mean. The spatial pattern of the averaged chill at each mode of PNA is shown in the Supplementary  Fig. 4, while the statistically significant differences in Fig. 7    the paired t-test are shown in Supplementary Fig. 5. The results show the PNA positive and negative modes have a similar spatial pattern of mode-averaged chill accumulation ( Supplementary Fig. 4). The spatial pattern of differenced results using CH (Fig. 7a-b) resembles those using CP (Fig. 7c-d). No matter which chill accumulation (either CH or CP) is considered, the focused region features a prevalent belowaverage CH (CP) at PNA positive mode (Fig. 7a or c), while a dominating above-average CH (CP) at PNA negative mode (Fig. 7b or d). Most of the differences during the PNA positive mode are not significant (Supplementary Fig. 5a and c), while during PNA negative mode, about 9% and 28% of the focused region shows significant chill accumulation above CH mean (Supplementary Fig. 5b) and CP mean ( Supplementary  Fig. 5d), respectively. Over the focused region, the significant CP difference during PNA negative mode had an average of about +2.1 CP above CP mean . Fig. 8 presents the difference between averaged CH and CP at each mode of PDO and CH mean and CP mean . The spatial pattern of the averaged chill at each mode of PDO is shown in the Supplementary Fig. 6, while the statistically significant differences in Fig. 8 that passed the paired t-test are shown in Supplementary Fig. 7. Again, there is no obvious difference of the mode-averaged chill accumulation between PDO positive and negative modes ( Supplementary Fig. 6). Fig. 8 indicates most of the focused regions present above-average chill accumulation for both CH and CP during the PDO positive mode and below-average chill accumulation during the PDO negative mode, though the results using CP ( Fig. 8a-b) are more prominent than those using CH (Fig. 8c-d). However, a significant difference is only observed when using CP (Supplementary Fig. 7c-d) rather than using CH ( Supplementary Fig. 7a-b). About 55% of the focused region presents significant above-average CP during PDO positive mode with an average of +2.3 CP above CP mean . Inversely, about 22% of the focused region exhibits below-average CP during PDO negative mode, with an average of −2.6 CP under CP mean ( Supplementary Fig. 7); this pattern is especially notable in the Sacramento Valley and south of the Central Coast.

Implication of decreasing trend of chill accumulation
Warming winter temperatures may have a significant impact on perennial crop production, particularly for specialty crops that require a necessary amount of chill accumulation to achieve optimum yields. The trend analysis of chill accumulation in this study (Fig. 3) shows that there is a statistically significant decreasing trend from 1979 to 2019 of chill accumulation (either using CH or CP) in the Central Valley and Central Coast of California. This result is consistent with previous studies which reported a declining winter chill in California from 1950 to the end of 21st century (Baldocchi and Wong, 2008;Luedeling et al., 2009). Based on observed temperature and climate projections, Baldocchi and Wong (2008) estimated the chill accumulation would be less than 500 chill hours per winter in the Central Valley of California by the end of the 21st century, while Luedeling et al. (2009) found that winter chill in the Central Valley may be decreased by 14%-33% by the same time. The research from Lobell and Field (2011) revealed that no currently-planted specialty crops will benefit from winter warming in California. In addition, Kerr et al. (2018) investigated the temperature sensitivity of selected specialty crops in California and found the specialty crops in San Joaquin and Sacramento Valleys present high winter sensitivity, which are vulnerable to a loss of winter chill. Although lowchill crops may continue to thrive in a warmer winter environment, adaptation efforts may be necessary for California to continue cultivating crops with high chill requirements.
In addition to increasing vulnerability and potentially shifting the geography of suitable cultivation locations (e.g., Parker and Abatzoglou, 2018), warming winter temperatures may impact the phenology and development of perennial specialty crops. Warmer winter temperatures have been shown to delay budburst and flowering in both natural (Yu et al., 2010) and agricultural systems (Asse et al., 2018;Luedeling and Gassner, 2012). Reduced chill accumulation may also result in poor fruit or nut development and reduced yield and quality (Atkinson et al., 2013;Luedeling et al., 2009;Saure, 1985). On the other hand, there is uncertainty in determining the required chill amount for specialty crops using different phenology stages. Pope et al. (2015) found the chill requirement of specialty crops at budbreak phase may be different from (higher than) the chill accumulation associated with greater potential yield, especially for pistachios and walnuts, meaning that the budbreak-based chill requirements may not reflect yield decline chill thresholds. Campoy et al. (2019) found the use of traditional approaches to determining chill requirements at 50% bud burst can result in an underestimation of the chill requirements to achieve optimal yield for tree crops (e.g. sweet cherry). The underestimation of chill requirements using traditional thresholds may aggregate the climate risk associated with warmer winters and subsequent declines in chill accumulation.

Implication of teleconnection-chill relationship
Although the impacts of large-scale teleconnection on California's climate, especially precipitation in southern California, have been widely studied (Schonher and Nicholson, 1989;Fierro, 2014;Allen and Luptowitz, 2017), less attention has been paid to integrated agriculture specific climate indices, such as chill accumulation. Jagannathan (2019) found that integrated climate metrics are more useful to farmers than basic climatic metrics for agriculture management. Jagannathan et al. (2020) also revealed that good skill in predicting basic climate metrics in California does not guarantee skill in prediction of decision-relevant metrics such as chill hours. Therefore, our work on assessing the impact of teleconnections on chill accumulation in California is necessary and important to support local agricultural decision-making.
This study makes two attempts to address the relationship between chill accumulation and teleconnection indices. First, the correlation analysis is conducted between monthly teleconnection indices and chill accumulation. The results reveal that PDO is positively correlated with CP accumulation with the strongest correlation (found in January), PNA presents a widely spread negative correlation with CP (with the most intense signals found in September), while ONI shows the weakest correlation with chill accumulation. This result is consistent with previous findings that PDO has a stronger impact on climate variability in wine regions in the western USA than ENSO (Jones and Goodrich, 2008); these wine regions include Napa Valley, which is located in our Central Coast focus region. Our results also showed that the correlation between chill accumulation and teleconnection indices is strongest during fall and winter, which agrees with previous findings that the signals of PDO and PNA are weaker during the summer than in other seasons (Mantua, 2001;van den Dool et al., 2000).
The second attempt is to compare the aggregated chill accumulation at different teleconnection modes (positive, negative, or neutral) with the long-term climatological mean of chill accumulation. The aim of this analysis is to tell farmers what to expect (above-or belowaverage chill) during different teleconnection modes. Our results reveal that teleconnection indices have a stronger response to CP variations than CH variations in California. The CP variation, especially the above-average signal, is most sensitive to the different modes of ONI, followed by those of PDO and PNA. For example, significantly aboveaverage CP can be expected during ONI positive mode (covering 57% area with an average of +3.6), PDO positive mode (covering 55% area with an average of +2.3), and PNA negative mode (covering 2.1% area with an average of +2.1). While significantly below-average CP is only observed during PDO negative mode, with an average of −2.6 covering 22% of the focused region.
Considering that chill accumulation is calculated using daily temperature, the impact of teleconnections on temperature variations can be used to explain the chill-teleconnection relationship. For example, Yu et al. (2012) revealed that the conventional Eastern-Pacific El Niño (corresponding to a positive mode of the ONI) is associated with negative winter temperature anomalies that are concentrated mostly over the southwestern states. This explains the significant above-average chill accumulation during ONI positive mode, considering that colder-thannormal temperature will result in more chill accumulation during winter.
Previous studies have also shown a coupling effect between ENSO and PDO. Sheppard et al. (2002) found ENSO and PDO can amplify each other resulting in increased climate variability over the US Southwest. Yu and Zwiers (2007) found the anomalous signals of climate change are stronger for the in-phase combinations of ENSO and PDO variability (El Nino/positive PDO and La Nina/negative PDO) than those for the out-of-phase combinations over North America. Birk et al. (2010) also revealed the ENSO-related variability was much higher during the PDO positive periods than those during the negative PDO periods. This may partially explain the above-normal chill accumulation observed during the positive modes of both ONI and PDO.
For PNA, previous studies have shown that the positive phase of the PNA is associated with above-average temperatures over western Canada and the extreme western United States (NOAA/National Weather Service, 2012). Abatzoglou and Redmond (2007) revealed that when the PNA has a negative trend in autumn this projects a cooling signal over western North America, and Abatzoglou (2011) found that positive values of the PNA are generally associated with anomalously high freezing levels, higher temperature and increased snowmelt across the Western United States. Further, Yu et al. (2012) found the impact of the Central-Pacific El Niño on US winter temperatures is associated with the PNA patterns, leading to a warmnorthwest, cold-southeast pattern of temperature anomalies. The positive correlation between PNA and temperature justifies the negative relationship between PNA and chill accumulation. All above findings confirm our results that PNA presents a significant negative correlation with chill accumulation, especially in September, and more chill accumulation is expected in negative PNA mode.
A better understanding of the relationship between teleconnections and chill accumulation could potentially help farmers in California to make careful considerations for future crop management. In the southeastern USA, the linkage between teleconnections and crop potential yield has been incorporated into a climate forecast information system (e.g., AgClimate) to assist agricultural decision-making (Fraisse et al., 2006). Similarly, the results and analysis from this study can be incorporated into and developed as an agricultural decision support tool in California, which would enable growers to manage risks under adverse chill conditions or take advantage of potentially favorable conditions.

Adaptation strategy
To cope with the decreasing winter chill, some management and adaptation strategies can be applied. Low chill years have been experienced in the past in California, especially in the southern San Joaquin Valley. In these cases, growers have applied rest breaking agents to reduce the negative effects of insufficient winter chill and encourage an even bloom (Glozer, 2006). Hussain et al. (2018) suggested seed treatments, application of plant growth regulators and compatible solutes, and use of plant mineral nutrients to cope with the insufficient chill. For existing orchards or vineyards, kaolin clay masks may be applied to increase chill accumulation, which have been found to increase chill accumulation by 5-7 CP during sunny, warm winters (Doll et al., 2015). Other possible strategies include genetic modifications, rootstock selection and developing cultivars with lower chill needs (Parker et al., 2020). For example, the southern highbush "low-chill" cultivars of blueberries only require only 150 to 600 chill-hours (Draper, 2007), making them promising cultivars for the San Joaquin Valley's mild winters (600 to 1200 chill-hours annually) (Bremer et al., 2008). Finally, changing crop management practices, such as moving crops to high elevations (>1000 m) to chill during dormancy, soil drying and pruning, may also help to strengthen the resilience to decreasing chill (Atkinson et al., 2013).

Broader impacts
The importance of sustainable agricultural development as a means of meeting global food and nutrition demands has never been more critical than in an era of climate change. The simultaneous consideration of socioeconomic, environmental and governance factors will be required. There is a recent policy shift from agricultural sustainability and food security via boosting agricultural production to a broader approach incorporating multidimensional targets, including nutritional effectiveness, social inclusion, and climate justice (Lang and Barling, 2012;Garnett et al., 2013;Agovino et al., 2018). The adjusted food sustainability proposed by Agovino et al. (2018) is a means to evaluate the policy efficiency on food sustainability of different countries. Policies are most effective when targeting nutritional challenges, while to some extent when focusing on agriculture. Other studies have shown (see Rusciano et al., 2020) the multidimensionality of sustainable development and indicate that a broader social-environmental innovation is promising to guide social, ecological and sustainable transformations.
Finally, sustainable development is a fundamental basis for climate justice, which can be conceived as the intertemporal climate equity and equality exchange among generations (Cisco and Gatto, 2021). Climate justice is closely connected with environmental, energy and resource justice (Schlosberg and Collins, 2014;Jenkins, 2018;Bickerstaff et al., 2013), and associated with reducing emissions, ensuring a just transition to a green economy, protecting vulnerable communities, and assuring intergenerational justice (Meyer and Roser, 2010;Shue, 2014). For example, in a recent study by Cisco and Gatto (2021) young generations may have to bear the brunt of sustainable development in the business-as-usual climate scenario and make more effort to reduce their leisure and consumption to reduce emissions, illustrative of the potential inequalities at play.
The rate of warming and other climate-mediated stressors challenges the sustainability and production demands for California agriculture. The reduction of chill requirements necessary for high value perennial commodities is just one, albeit important, aspect of climate change. Agricultural operations will be challenged by the increased frequency of extreme heat, water scarcity, and pest and pathogen pressure. Ameliorating the future impacts of climate change will require socio-economic and technological innovations, and science-based information to support climate-informed policy and decision making at local and regional scales (e.g., counties, water districts, state government). Such an approach could yield robust sustainability indices with broad implications including those related to disadvantaged communities, food and nutrition security, and environmental justice. After all, humans are the foundation of sustainable agroecological systems and our ability to adapt, in the most comprehensive sense, needs to be part of the conversation.

Limitations and future work
The findings from this study should be placed in context by acknowledging the study limitations, including: (1) Only one month of teleconnection indices is used to classify different modes. Although the month with the maximum correlation is used, this is still a simplified approach and the results may vary considering the lagged correlation. A more robust method is to use the average of teleconnection indices over multiple months, which can be explored in future studies.
(2) The PDO requires a long data record to compute (e.g., more than 50 years). However, the data period in this study is limited to 41 years (1979 to 2019). The results on PDO may be more robust when longer time series data are used.
Nevertheless, this study laid a solid foundation for our next step of dynamic forecasting of chill accumulation. The teleconnection indices at the months where maximum correlations are achieved can be used as the predictor to forecast the chill accumulation well before the chill accumulation period begins. Statistical models (e.g., multilinear regression or stepwise regression), as well as machine learning methods (e.g., random forest, gradient boosting), can be used and compared in the prediction of chill accumulation. Focus can be given to specialty crops, and crop-specific models can be developed. The early release of forecasted chill accumulation can be helpful for growers in planning their cropping and management activities well before the chill season begins.

Conclusions
Based on the above analysis, we came to the following conclusions.
(1) The trend analysis shows that there is a statistically significant decreasing trend from 1979 to 2019 of chill accumulation (either using CH or CP) in the Central Valley and Central Coast of California. The decreasing rate is −4.98CH/year and −0.15CP/ year, respectively. Multiple management and adaptation strategies have been introduced to cope with the decreasing winter chill.
(2) In general, CP presents a stronger response to teleconnection indices than CH in California. In terms of correlation, ONI has the weakest impact among the three teleconnection indices and explains on average of 11% CP variation. PDO has the strongest correlation among the three, it explains on average of 16% CP variation, ranging between 10% to 39%. PNA shows the most widely spread negative correlation with CP, covering 85% of the focused region and explaining on average of 15% CP variation. The correlation signals are weaker during the summer than other seasons, and the strongest correlation can be observed well before the start of the chill season. The dominant months when the maximum correlation is achieved using either CH or CP are September for ONI, January for PDO and September for PNA. This provides a potential opportunity for early-season chill prediction using teleconnection indices. (3) When aggregated according to different teleconnection modes, the CP variation, especially the above-average signal, is most sensitive to the positive modes of ONI, followed by the positive mode of PDO and negative mode of PNA. During ONI positive mode, +3.6 above-average CP can be expected; +2.3 above-average CP is expected during PDO positive mode, while +2.1 above-average CP can be expected during PNA negative mode. The results of this analysis can tell farmers what to expect (above-or belowaverage chill) during different teleconnection modes, and make them better prepared for crop management. (4) Agricultural management and adaptation may benefit from the direct investigation of teleconnection impacts on integrated climate variables (e.g., chill accumulation). Better understanding the relationships between teleconnections and chill accumulation could potentially help farmers in California to make careful considerations for future crop management. The methods from this study can be adapted for other regions, and the analysis from this study can be easily incorporated into and developed as an agricultural decision support tool to help growers and policymakers to make adaptive management decisions.

Funding sources
This work was supported by the funding from USDA-ARS and The USDA California Climate Hub, Award number: 58-2032-8-059 and University of California Office of the President, Multicampus Research Programs and Initiatives funding.