Spatiotemporal distribution of the potential risk of frost damage in tea fields from 1981 - 2020: A modeling approach considering phenology and meteorology

Climate change may induce severe frost damage to crops, and thus a reasonable assessment of frost risk, considering both crop phenology and meteorology, is required. Here, we introduced a new index of potential frost risk ( F risk ) using thermal time ( minimum air temperature below the threshold value ) weighted by the percentage of budburst ( P bud ) . Moreover, we evaluated the spatiotemporal distributions of F risk in tea fields within a 60 km × 60 km area in east Japan from 1981 - 2020, using 1 km 2 -gridded meteorological data and a newly developed model of P bud . The P bud model considered three phenological phases ( endodormancy, ecodormancy, and progress of budburst ) and successfully represented changes in the P bud of the tea buds for 15 years, with root mean square errors of 8.5 percentage points. The spatiotemporal distributions of F risk over the past 40 years showed that potential frost risk significantly increased at elevations ranging from 50 m to 300 m because the budburst advanced at a faster rate than the temperature warming. These elevations corresponded to areas where tea plants were mainly cultivated, which indicates that tea cultivation is becoming vulnerable to frost, and the risk of economic losses due to the frost is increasing. The proposed assessment of frost risk could contribute to predicting frost damage and developing more reliable strategies for the operation of frost protection under the effects of future climate change.


Introduction
Frost often damages crops and induces huge economic losses in agriculture Snyder et al., 2005 . In recent years, tremendous frost damage has been reported in fruit orchards in the United States 112 million USD in 2007; Gu et al., 2008 and in tea fields in China 268 million USD in 2010; China Tea Marketing Association, 2010 and Japan 45 million USD in 2010; Matsuo et al., 2010 . Owing to recent climate warming, the risk of such frost damage may be expected to decrease in the future because of advances in the date of the last spring frost Schwartz et al., 2006 . However, climate warming has also caused advances in phenological events such as budburst, leafout, and flowering in spring Parmesan and Yohe, 2003 , reducing the cold hardiness of plants. These advances in spring phenology and extreme temperature events, such as sudden cold spells in future climates, may increase the risk of frost damage Augspurger, 2013;Ma et al., 2019 . Many investigations into frost risk have been performed based on trends in observed minimum temperatures, which give the dates of freezing events and the resulting length of the growing season. By analyzing data of daily minimum temperatures observed at ~1,000 meteorological stations, Easterling 2002 found that the number of frost days decreased in the United States from 1948-1999. Robeson 2002 used data from a longer observation of minimum temperature in Illinois from 1906-1997 and found trends of the earlier spring freezing events. Some investigations have considered plant phenology as well as minimum temperature trends. To quantify the potential frost risk, a false spring index FSI; Marino et al., 2011 is defined as the difference between the date of the last spring freeze and the date of phenological events such as budburst, leafout, and flowering, and positive values of FSI indicate risk situations. Vitasse et al. 2018 evaluated the FSI of tree species in Switzerland and revealed that the frost risk increased at higher elevations and remained the same at lower elevations from 1975-2016. The use of FSI is a reasonable approach that accounts for both meteorological and phenological aspects. However, FSI depends on the date when a certain percentage of a population reaches the stage of budburst, leafout, and flowering, and thus ignores temporal changes in the percentage of phenological events in the population. As the percentage increases with time i.e., the number of buds, leaves, and flowers that have lost cold hardiness increases , the potential frost risk and resulting agricultural losses probably increase. Thus, considering the changes in the percentage of phenological events can lead to a more reliable assessment of frost risk.
Phenological events are often estimated using an empirical model because of limited data on observed phenology, and many models have been developed to estimate phenological development in crops e.g., Horie and Nakagawa, 1990, Nakagawa and Horie, 1995and Yin et al., 1997Weir et al., 1984, Maruyama et al., 2010and Kawakita et al., 2020Williams et al., 1985, Caffarra and Eccel, 2010and Sugiura et al., 2019 in grapes . Conventional phenological models of budburst, leafout, and flowering, which are required for the assessment of frost risk, are mainly divided into two types: one-phase models e.g., Cannell and Smith, 1983;Chuine and Cour, 1999 , considering the ecodormancy processes affected by forcing temperatures in spring, and two-phase models e.g., Richardson et al., 1974;Chuine, 2000 , which consider endodormancy processes affected by chilling temperatures from autumn to winter as well as the ecodormancy processes. These models generally predict the dates of phenological events not the percentages , and require some modification to evaluate the changes in the percentage of phenological events.
Tea plants Camellia sinensis L. O. Kuntze , one of the most important economic crops around the world, have been seriously damaged by recent frost, leading to huge economic losses as mentioned above. Given the critical importance of assessing frost risk under current warming climate trends, some researchers have considered minimum temperatures in their investigation of the frost risk in tea fields e.g., Li et al., 2018 , and a model to predict phenology of tea buds cold acclimation and deacclimation has been developed Kimura et al., 2021 . However, frost risk in terms of both meteorology and phenology, remains unclear in tea plants. Tea plants are cultivated in a wide range of regions, from plains to mountains, and frost risks vary in space and time at different scales in different species Bennie et al., 2010;Vitasse et al., 2018 . Thus, the spatial variability of frost risk should be assessed in terms of reasonable strategies for the application of frost protection e.g., frost protective fan, Kimura et al., 2017;and sprinkler, Suzuki et al., 1993 across a wide range of regions and for a long period under future warming climate.
The objectives of this study were 1 to develop a new index for the potential frost risk based on changes in the percentage of phenological events and in daily minimum temperature, 2 to model the changes in percentage of the phenological event in this study, budburst in tea plants , and 3 to assess the spatiotemporal distribution of the potential frost risk in tea fields using the developed index and the model over the period 1981-2020 on a regional scale.

Study area
A 60 km 60 km area located in east Japan from 35.78 to 36.25 N latitude and 138.92 to 139.63 E longitude was selected Fig. 1a, b to assess the spatiotemporal distribution of the potential frost risk. In this area, elevation ranges from 0 to 1,300 m Fig. 1b , and tea plants are widely cultivated from plains to mountains Omata, 2008 . Representative points P1 -P4 in Fig. 1b near automated meteorological data acquisition systems AMEDAS were selected to compare long-term trends in the potential frost risk at different elevations. The study area is the northern limit of economical tea production in Japan and has a humid subtropical climate. According to averaged records of AMEDAS stations within the area 17 stations from 1981-2019, an annual mean temperature of 14.2 C and an annual mean precipitation of 1,427 mm have been recorded in the area.

Evaluating potential frost risk
We developed a new index for the potential frost risk during the growing season F risk under the assumption that a higher percentage of a phenological event in a population results in a larger risk of frost damage under the same temperature. This assumption was based on the finding that cold hardiness tended to decrease gradually during bud development in tea plants Yanase, 1972 and other tree species Lenz et al., 2013 . Using this assumption, F risk can be defined as the accumulation of frost risk at a given date t F risk,t , which is expressed by thermal time the daily minimum air temperature [T min,t ] below the threshold value [T frost ] weighted by the percentage of budburst P bud,t . The F risk,t shows positive values i.e., risk situations only if T min,t is below T frost and P bud,t is above 0 i.e., budburst conditions . The calculation of Eq. 1 began when budburst definitely did not occur D s = 1 in day of year within the study area and continued until T min was definitely not below T frost D e = 150 in day of year within the study area. In this study, T frost was set at 3 C, where frost protection e.g., frost protective fans is conventionally required in tea fields Araki et al., 2010 .

Modeling percentage of budburst
For modeling P bud , the phenological cycle for tea budburst is divided into three phases: endodormancy from autumn to winter, ecodormancy from winter to spring, and initial budburst and subsequent increase of P bud Fig. 2 .
Endodormancy, which coincides with the cold acclimation of buds, was assumed to be induced by decreasing temperature and released by fulfilling the chilling requirement or increasing the photoperiod Charrier et al., 2011;Laube et al., 2014;Du et al., 2019 . Chilling requirement at a given date t CR t was based on the cold acclimation model imitating cold stress memory in plants Kimura et al., 2021 , which was expressed as accumulation of thermal time daily average air temperature [T ave ] below the threshold value [T ac ] corresponding to the temperature when endodormancy starts in autumn weighted by a forgetting curve: where k and c are the parameters of the forgetting curve, and endo is the date when endodormancy starts. The forgetting curve shows values from 0 to 1, and the values are highest on the latest day, decreasing gradually as the number of days decreases Fig. 2 . Thus, Eq. 3 shows that the thermal time in the latter days n close to t strongly contributes to cold acclimation and that in the former days n close to endo , it does not Kimura et al., 2021 . The calculation of Eq. 3 began when T ave first fell below T ac in autumn within the study area. When CR t reached the critical value CR crit , endodormancy was released and the ecodormancy phase started. If the chilling requirement did not reach CR crit , the critical daylength DL crit offset insufficient chilling, and the ecodormancy phase started Kimura et al., 2021 . Ecodormancy, which coincides with cold deacclimation of buds, was assumed to be driven by competition between chilling and forcing requirements Chuine, 2000 at a given date t CF t as modeled by Kimura et al., 2021 : where FR t is the forcing requirement at a given date t, which was calculated by the thermal time above the threshold temperature T de , and eco is the date when ecodormancy starts. The calculation of Eq. 5 began when the ecodormancy phase started. Ecodormancy was released when CF t reached the critical value CF crit , and when the initial budburst occurred i.e., P bud was first above 0 . After the initial budburst, P bud increased in response to warm temperatures, so that the daily change in P bud was modeled by:

Fig. 2.
Modeling of the percentage of budburst P bud , divided into three phases: endodormancy, ecodormancy, and increase in P bud . Details of equations and model parameters are described in section 2.3 and Table 1.
where a is the parameter of the temperature response curve for P bud , and T 50 is the temperature where P bud increased by 50 percentage points pp. in a day. The calculation of Eq. 6 continued until P bud reached 100 and P bud remained constant at 100 .

Model parameterization and validation
For parameterization and validation of the model, P bud was observed in a tea field located in Saitama, Japan Saitama Tea Research Institute; 35 4829 N, 139 2052 E, 148 m a.m.s.l., Fig. 1b . The observation was conducted for the cultivar Yabukita, which is widely cultivated in Japan, from mid-April to late-May for 15 years 2006-2020 , at intervals of ~5 d. Budburst was defined as when the length of the new bud enveloped by leaves was twice as long as these leaves The annual scientific conference of tea, 1986 , and P bud was calculated as the ratio between the number of budbursts and the number of non-bursting buds n = ~100 leaves .
The model parameters were assessed by minimization of the root mean square error RMSE between the estimated and observed values of P bud for 15 years. For the parameterization, the differential evolution algorithm was applied using the 'DEoptim' R package version 2.2-5 , as proposed by Mullen et al. 2011 . A two-fold cross validation was conducted. The dataset for 15 years was randomly partitioned into two sets d 1 for 8 years and d 2 for 7 years , then the model was parameterized on d 1 and validated on d 2 , followed by parameterizing on d 2 and validating on d 1 .
To reduce model complexity, three parameters T ac , T de , and c were fixed at the values established by Kimura et al. 2021 from 10 years of observed cold-hardiness data from the same cultivar Yabukita. These fixed parameters were determined on the basis of the Akaike information criterion AIC; Akaike, 1974 , which can consider both model accuracy and complexity. Among various combinations of fixed and fitted parameters, the combination was chosen wherein AIC showed the smallest value i.e., the combination in which T ac , T de , and c were fixed, and other parameters were fitted . The optimized parameters were determined using all data for 15 years Table 1 .
To calculate P bud and F risk , meteorological data of T min and T ave during 1981-2020 within the study area were obtained from the Agrometeorological Grid Square Data AMGSD; Ohno et al., 2016 . The AMGSD provides 1 km 2 -gridded data of T min and T ave using spatial interpolation among the observed values at AMEDAS located at intervals of ~20 km in Japan and corrects for the differences in temperature due to elevation. These data were used to study the effects of climate change on crops e.g., Nakano et al., 2017;Fukui et al., 2019 . Photoperiods were calculated using the day of the year, sun position, and latitude in each grid.

Spatiotemporal analysis of potential frost risk
The values of P bud , T min accumulated T min below T frost [ 3 C ] after initial budburst , and F risk from 1981-2020 were assessed in each grid of AMGSD within the study area. The significance of trends in P bud , T min , and F risk from 1981-2020 were determined using non-parametric Mann-Kendall tests at 5 significance level Kendall, 1975 , and then Sen's slope Sen, 1968 , which represents the degree of trends, was assessed in the grids wherein significant trends were found. The grids wherein T min was below 15 C in winter were removed from the trend analysis because the maximum cold tolerance in Yabukita is around 15 C in midwinter Kimura et al., 2021 .

Model validation
Seasonal changes in T ave , which is the input variable of the model of P bud , showed different patterns over 15 years Fig. 3 . In 2011-2012, T ave represented the coldest winter and spring, and, in contrast, 2019-2020 represented the warmest winter and spring. In 2012-2013 and 2017-2018, the changes in T ave showed an alternative pattern, that is, colder autumn-winter and warmer spring.
These different changes in T ave resulted in different durations of endodormancy and ecodormancy Fig. 3 and large variations of P bud Fig. 4 in the model simulation. Except for the initial budburst in 2007, the present model of P bud successfully captured the observed variations of P bud for 15 years Fig. 4 , even when tea buds were exposed to both colder and warmer climates. The error variance in the validation dataset was higher in the high P bud region. Overall, the observed variation of P bud explained 91 and 87 of the predicted variation of P bud for the calibration and validation datasets, respectively, and the accuracy of the model was less than 10 . root mean square error [ RMSE ] = 7.0 pp. in calibration and 8.5 pp. in validation Fig. 5 .

Spatiotemporal distributions of budburst, minimum
temperature, and potential frost risk Trends in P bud , T min , and F risk were divided into two types of patterns Fig. 6 . In the plain P1 and mountain P3 areas, Temperature where percentage of budburst increases by 50 pp. in a day 33.6 C Fitting this study the date of budburst advanced over the course of 40 years 0.18 d y -1 in P1 and 0.30 d y -1 in P3, for the initial budburst , but T min did not increase. Consequently, F risk did not increase in these areas. By contrast, in the hilly areas P2 , the date of budburst advanced 0.33 d y -1 for the initial budburst and T min increased 0.16 C y -1 , together leading to a significant increase in F risk . An exceptional trend was found in P4 Fig. 7 . Significant advances in P bud were not detected at P4, but T min significantly increased over the 40 years. This increase in T min resulted in a significant increase in F risk .
The dates of the initial budburst significantly advanced across almost all study areas Fig. 8 , and the areas where significant advances were found decreased with the increase in P bud . The advance rates were higher in the hilly and mountainous areas than in the plain areas and decreased with increasing P bud across the study area Fig. 8 .
Significant increases in T min were concentrated in the hilly areas Fig. 9 . In almost all these areas, the increase in T min was attributed to the advance of budburst, but in the limited area near P4, T min increased even though the dates of budburst did not advance. The spatial distribution of the trend in F risk approximated that of T min Fig. 8 , with slight differences because T min was biased by P bud for evaluating F risk . Overall, the significant increase in F risk was concentrated at elevations ranging from 50 m to 300 m within the study area Table 2 . These elevations corresponded to hilly areas where tea plants were mainly cultivated.

Evaluation of F risk and P bud model
To date, potential frost risk has often been evaluated by comparing dates between the last freeze and budburst Marino et al., 2011 , as mentioned earlier. However, to assess frost risk, the period between initial budburst and full leafout the period of vulnerability to frost should be considered as well as the date of budburst Chamberlain et al., 2019 . Thus, F risk was formulated in this study by introducing the thermal time T min below T frost weighted by P bud . To our knowledge, this study is the first to consider the impact of the changes in budburst percentage on frost risk. The F risk calculation assumed that higher P bud in a population has a larger frost risk; this assumption could be suitable because cold hardiness tends to decrease gradually during bud development in tea plants Yanase, 1972 and in the other tree species Lenz et al., 2013 . There have been many models for predicting the dates of phenological events e.g., Chuine, 2000;Caffarra et al., 2011 ; however, to our knowledge, there is no published model for predicting P bud with consideration to endodormancy and ecodormancy phases. Conventional models of budburst have High accuracy of the model for P bud was noted with the validation using 15 years of data, including unusually cold and warm weather, indicating that this model could predict the development of tea buds that experienced different temperature changes. The fixed parameter T ac had the value of 17.4 C, indicating that endodormancy cold acclimation in tea plants started when the daily average temperature fell below 17.4 C. At a similar temperature, an increase begins in tea plants of sugar levels and sugar-related gene expression, which play essential roles in driving cold acclimation Yue et al., 2015 . Eq. 6 in the P bud model the function newly introduced to express the increase of the budburst percentage produced a logistic curve of increasing P bud even in a high temperature region, although plant growth is generally suppressed during particularly high temperatures Warrington and Kanemasu, 1983 . Moreover, Eq. 6 had the fitted parameter T 50 of 33.6 C. This value indicates that 33.6 C in daily average temperature increases budburst by 50 pp. in a day, but validity of the fitted value was not investigated in the present study. The logistic function and the fitted parameter would be applicable within the study area, but a large variance across high-P bud regions in the model validation may be attributed to the function and the parameter value. Thus, Eq. 6 should be validated in future studies by measuring P bud in different locations with a wide range of temperature conditions.

Spatiotemporal distributions of frost risk in the tea fields
A few investigations of frost risk have been conducted in tea fields. Li et al. 2018 evaluated the risk of tea frost in Zhejiang in China from 1971 to 2017, using daily minimum temperature data, a digital elevation model, and slope data. They indicated that the risk was higher in hilly and mountainous areas than in plain areas, finding similar results to the present study. Lou et al. 2013 evaluated the risk of tea frost in Longjing in China from 1971 to 2010, as the difference between the beginning date of tea plucking and the date of last frost. They indicated that almost all the points 12 meteorological stations had no significant trends in the risk, but this result cannot be directly compared to our results because data on the elevation was not provided in their study.
This study consistently highlighted the significant increase in frost risk in hilly areas 50 m < elevation < 300 m , where tea plants are mainly cultivated. This is because the budburst advanced at a faster rate than the temperature warming in the areas, and, consequently, the period of budburst further overlapped the period when T min was below T frost 3 C see T min at P2 in Fig. 6 . These results indicate that tea plants in the main cultivation areas are becoming vulnerable to frost, and the risk of economic losses is likely to increase. By contrast, a trend in the frost risk in the mountainous areas elevation > 450 m was not found because the date of the initial budburst rarely reached the period where T min was below T frost 3 C over the period of 1981-2020 see T min at P3 in Fig. 6 , despite the faster advancement of budburst in the mountainous areas than in the hilly areas. The higher rate of advancement of budburst at higher elevations was attributed to the amplification of the warming rate with elevation Pepin et al., 2015 . This elevation-dependent warming EDW impacts changes in ecosystems and farming communities in mountainous environments Pepin et al., 2015 . In this study, F risk in the mountainous areas did not increase over 40 years, but if EDW continues, F risk may increase rapidly in the near future. The increased risk of frost exposure at high elevations has been demonstrated in tree species in Switzerland from 1975-2016Vitasse et al., 2018 In the plain areas elevation < 20 m , no significant trend in frost risk was attributed to the brief vulnerability to frost, as is the case for the mountain areas.

Limitations of this study
One significant limitation of this study is that F risk does not consider the relative contributions of T min and P bud contributed to the potential frost risk; this is because F risk is formulated by simple multiplication Eq. 2 . For a more reliable assessment of frost risk, the contribution of T min and P bud should be explicitly quantified by measuring the actual frost damage. In tea fields, frost protective fans are usually operated to prevent frost damage Kimura et al, 2017 . The frost protective fan mitigates temperature inversion on radiative frost nights and increases the surface temperature of tea plants, probably leading to some temperature difference between a meteorological station and an actual tea field. For more accurate evaluation of the spatiotemporal distributions of F risk and P bud , the temperatures in the actual tea field should be used. However, this is hampered by the necessity for sensor networks consisted of numerous accurately corrected thermometers.
The applicability of the bud model to different cultivars and other climate zones was unclear. Particularly, CR crit and CF crit , which are temperature requirements determining the transition of the phenological phases in the model, may vary with latitude even in the same species McKown et al., 2018 , under the influence of such other factors as precipitation and the resultant changes in soil moisture. To develop a more robust model and apply the model to wider regions, a set of data on the P bud observed for different cultivars and in other climate zones is required.
The exceptional trend near P4 indicates that the period of vulnerability to frost increased despite the absence of a significant advance of budburst; that is, the minimum temperature decreased over the past few decades. This may perhaps be caused by changes in the local geographical effects of urbanization; land use near P4 has tended to urbanize over the past 15 years, according to the National Land Numerical Information provided by the Japanese Ministry of Land, Infrastructure, Transport, and Tourism.
Finer spatial resolutions of temperature are needed to more Fig. 7. Trends in day of year DOY in which percentage of budburst P bud exceeded 0 initial budburst and 50 and reached 100 , accumulated daily minimum temperature below 3 C after the initial budburst T min , and potential frost risk F risk from 1981-2020 at the representative points P4. The location of P4 is described in Fig. 1. Linear lines represent Sen's slopes with Mann-Kendall tests at 5 significance level, and data without lines show no trends at 5 significance level. and reached 100 , accumulated daily minimum temperature below 3 C after the initial budburst T min , and potential frost risk F risk from 1981-2020 at the representative points P1, P2, and P3. The locations of P1 P3 are described in Fig. 1. Linear lines represent Sen's slopes with Mann-Kendall tests at 5 significance level, and data without lines show no trends at 5 significance level. meaningfully assess distributions of F risk and P bud . In the hilly and mountainous areas where the tea plants are usually cultivated, significant differences in temperature occur within a 1 km 2 grid section, mainly due to cold air drainage on radiative frost nights Chung et al., 2006 , leading to a bias of P bud within the grid. Chung et al. 2006 developed a minimum-temperature prediction model which considered cold air drainage and temperature inversion strength in complex terrain, using a digital elevation model. Moreover, Maclean et al. 2019 provided a method and an R program to assess the spatially fine distribution of air temperature resolution < 30 m , considering the local effects of cold air drainage and other microclimate conditions. These methods were   . 9. Spatial distributions of Sen's slope degree of trend of accumulated daily minimum temperature below 3 C after the initial budburst T min and potential frost risk F risk , from 1981-2020. Blank areas show no trends at 5 significance level using Mann-Kendall tests. The representative points P1 P4 are indicated on the map. based on the calculation of airflow direction and accumulation within a study area, allowing the assessment of where cold air converges. The integration of these methods with the bud model could provide more meaningful assessments of frost risk.
Nevertheless, we believe that the new index of potential frost risk and the model of budburst percentage proposed in this study could improve the conventional assessment of frost risk and contribute to the development of reliable strategies for the application of frost protection under the effects of future climate change.

Concluding Remarks
A new index for potential frost risk, considering phenology and meteorology, was proposed using thermal time minimum air temperature below the threshold value weighted by the proportion of budburst. The proportion of budburst in tea plants was modeled, by formulating endodormancy, ecodormancy, and the subsequent increase of the number of budbursts. This model successfully predicted the proportion of budburst in tea plants over 15 years, with root mean square errors of 8.5 percentage points. Using the proposed index and the model, spatiotemporal distribution of frost risk in tea plants was evaluated in east Japan from 1981-2020. The distribution showed that frost risk increased over that period in hilly areas where tea plants are mainly cultivated. This result suggests an urgent need for improved frost protection in tea fields.