Mapping Threats of Spring Frost Damage to Tea Plants Using Satellite-Based Minimum Temperature Estimation in China

: Spring frost damage (SFD), deﬁned as the disaster during the period of newly formed tea buds in spring caused by lower temperature and frost damage, is a particular challenge for tea plants ( Camellia sinensis ), whose capacity to adapt to extreme weather and climate impacts is limited. In this paper, the region of the Middle and Lower Reaches of the Yangtze River (MLRYR) in China was selected as the major tea plantation study area, and the study period was focused on the concentrated occurrence of SFD, i.e., from March to April. By employing the standard lapse rate of air temperature with elevation, a minimum temperature (Tmin) estimation model that had been previously established was used based on reconstructed MYD11A1 nighttime LST values for 3 × 3 pixel windows and digital elevation model data. Combined with satellite-based Tmin estimates and ground-based Tmin observations, the spatiotemporal characteristics of SFD for tea plants were systematically analyzed from 2003 to 2020 in the MLRYR. The SFD risks at three scales (temporal, spatial, and terrain) were then evaluated for tea plants over the MLRYR. The results show that both SFD days at the annual scale and SFD areas at the daily scale exhibited a decreasing trend at a rate of 2.7 days/decade and 2.45 × 10 4 ha/day, respectively (signiﬁcant rates at the 0.05 and 0.01 levels, respectively). The period with the highest SFD risk appeared mainly in the ﬁrst twenty days of March. However, more attention should be given to the mid-to-late April time period due to the occurrence of late SFD from time to time. Spatially, areas with relatively higher SFD days and SFD risks were predominantly concentrated in the higher altitude areas of northwestern parts of MLRYR for both multi-year averages and individual years. Fortunately, in regions with a higher risk of SFD, the distribution of tea plants was relatively scattered and the area was small. These ﬁndings will provide helpful guidance for all kinds of people, including government agencies, agricultural insurance agencies, and tea farmers, in order that reasonable and effective strategies to reduce losses caused by spring frost damage to tea plants may be recommended and implemented.


Introduction
In early spring, herbaceous and woody plant species begin to sprout. It is well known that the meteorological conditions during the flowering and fruit setting periods for economic fruit trees are critical for determining the output and economic income. Similarly, the meteorological conditions are also important for tea (Camellia sinensis) plant production during both the newly formed tea bud period and the spring tea picking period, because the economic benefits of spring tea have accounted for more than 70% of the total tea revenue in recent years in the middle and lower reaches of the Yangtze River region (MLRYR), China [1]. As one of the major tea plantation areas in China, the MLRYR tea plantation methods [24,25]. These methods have been proven to be successful in estimating temperatures near meteorological stations. However, the spatial variability of temperature may be sensitive to latitude, longitude, elevation, and season. A useful interpolation method should be able to capture characteristics in different research regions and periods by adding appropriate variables [26]. Unfortunately, such a universal interpolation method has not been effectively developed yet.
Satellite remote sensing provides unprecedented globally continuous coverage of critical land surface parameters at high spatial resolution. These parameters include vegetation indices, leaf area index, surface temperature, etc. Land surface temperature (LST) is the radiative skin temperature of the land surface [27] as retrieved from top-of-atmosphere brightness temperatures in the infrared spectral channels of a remote sensor. Spatial estimates of temperature at high temporal (daily) and spatial resolution (1 km) appeared in 1981 with the launch of the Advanced Very High-Resolution Radiometer (AVHRR) on board a National Oceanic and Atmospheric Administration (NOAA) satellite, and later with the Moderate Resolution Imaging Spectroradiometer (MODIS) on board NOAA's Aqua and Terra satellites in 2000 [28]. Therefore, methods regarding air temperature (Ta) estimation with spatial patterns have been widely researched based on MODIS LST products and some auxiliary data over wide areas, such as the mixed bootstrap and jackknife resampling method in Portugal [29], the stepwise linear regression method in east Africa [30], the temperature-vegetation index (TVX) method in the north Tibetan Plateau [21], the random forest machine learning approach in urban landscapes [31].
Methods for estimating maximum (Tmax), minimum (Tmin), and average (Tavg) air temperature have been successfully developed in recent years. However, applications of spatially continuous Ta patterns from remote sensing data have rarely been reported in publications. Wang et al. took three cooling periods in 2006 as examples to identify the spatial distribution of SFD for tea plants with estimated Tmin based on MODIS LST products and the standard lapse rate of air temperature with elevation [32]. Results showed that three cooling periods were exactly identified, especially for several highlighted cities and counties recorded in the "Yearbook of Meteorological Disasters in China".
The major objective of this study was to further analyze the spatiotemporal distribution characteristics of SFD on tea plants with spatially continuous satellite-based Tmin estimates at high temporal scale over the MLRYR in China from 2003 to 2020 based on a previous publication [32], and then to clarify the temporal, spatial, and terrain risks of SFD. To reach this goal, the study area and data are first described in Section 2. The Tmin estimation method based on satellite-based reconstructed LST values and elevation products is then briefly introduced in Section 3. Next, the spatiotemporal distribution characteristics and risks of SFD for tea plants are systematically analyzed at 1 km resolution over the MLRYR in China in Section 4. Finally, some limitations are discussed and conclusions are summarized at the end of this paper.

Study Area
The MLRYR (Figure 1) is across six provinces (Jiangsu Province, Zhejiang Province, Anhui Province, Jiangxi Province, Hubei Province, Hunan Province) and one municipality (Shanghai Municipality). The MLRYR is the most important tea-production area in China, and its annual tea output accounts for about two thirds of the country's total output [12]. The tea gardens in the MLRYR are mainly distributed in hills and low mountain areas, and a few are scattered in high mountain areas with the elevation being more than 1000 m above sea level ( Figure 1). The climatic characteristics are seasonally distinguishable in the MLRYR, and are characterized as warm in spring, hot in summer, cool in autumn, and cold in winter. , the elevation in four levels was mapped with different shades of yellow, and the land cover type of shrub was denoted with red dots.

Ground-Based Meteorological Data
Ground-based daily Tmin observations at 163 meteorological stations from March to April (the period when SFD occurs for tea in the MLRYR) in 2003-2020 were downloaded from the China Meteorological Administration Meteorological Data Network (http://data. cma.cn/site/index.html. Last visited on 25 April 2021). All data were quality checked, and there were no missing temperature records during the study period in the MLRYR.

Remotely Sensed Products
(1) MODIS LST products MODIS products are available for public download via the Earth Observing System Data and Information System (EOSDIS) of the National Aeronautics and Space Administration (NASA (https://earthdata.nasa.gov/. Last visited on 26 April 2021). In this study, daily MODIS LST products (MYD11A1) at a spatial resolution of 1 km on board the Aqua satellite platform were used. To cover the entire study area, the LST products of four tiles in a 1200 by 1200 grid (h27v05, h27v06, h28v05, and h28v06) were downloaded from March to April in 2003-2020. A previous study demonstrated that nighttime MYD11A1 LST averaged at 3 × 3 pixel-window resolution was the best proxy for daily Tmin estimation [32]. Therefore, four tiles each day were re-projected to Albers Conical Equal Area (ACEA) projection, and then merged, resized, and mosaicked to get a nighttime LST dataset that was exactly consistent with the study area. Then daily good retrieval LST values recognized by quality control (QC) flags within a moving window of 3 × 3 pixels were averaged for estimating daily Tmin.
(2) MODIS Land Cover Type products The MODIS Land Cover Type Product (MCD12Q1) was created using a hierarchical classification model that ingested a full year of 8-day MODIS Nadir BRDF-Adjusted Reflectance (NBAR) data [33]. The MCD12Q1 product provides six different legacy classification schemes, including IGBP (International Geosphere-Biosphere Programme), UMD (University of Maryland), LAI (Leaf Area Index), BGC (BIOME Biogeochemical Cycles), PFT (Plant Functional Types), and LCCS (Land Cover Classification System), at annual time steps and 500 m spatial resolution for 2001-present [33]. In this study, the PFT scheme was used to extract the shrub classification that corresponded to the tea planting area in the MLRYR (red color dots in Figure 1b). Data can be freely downloaded at a 1 • × 1 • tile for convenience of distribution and handling by users. The total 168 tiles covering from 24 • to 35 • N latitude and from 109 • to 122 • E longitude were merged, resized, and mosaicked, and the elevation values within a moving window at 3 × 3 km grid cells were aggregated up to get a DEM image with the same spatial resolution of 1 km for all remotely sensed data.
All remotely sensed products were converted to *.tif format files by MODIS Reprojection Tools (MRT). Spatial calculation in pixel-by-pixel was done in Interactive Data Language (IDL), and spatial exhibition and mapping were drawn in ArcGIS. Additionally, the excel software was also used in statistical analysis.

Methodology
The threats of SFD to tea plants were mapped using satellite-based Tmin estimation in the MLRYR in China, and its pipeline was shown in Figure 2. The data used in this study includes ground-based meteorological observations (describes in Section 2.2.1) and remotely sensed products (in Section 2.2.2). The methodology mainly concerns three parts in this pipeline: (i) satellite-based Tmin estimation (in Section 3.1); (ii) spatiotemporal characteristics of SFD expressed as SFD days and SFD stations with ground-based Tmin observations, and SFD days and SFD pixels with satellite-based Tmin estimations (in Section 3.2); and (iii) three-dimensional risks of SFD with pixel-based SFD percentages and SFD frequencies (in Section 3.3). The spatiotemporal characteristics of SFD (in Sections 4.1 and 4.2) and its risk (in Section 4.3) to tea plants were exhibited in the bottom of the pipeline. Here, three gray boxes represent satellite-based data or method to distinguish them from the ground, and the blue number signs (such as Sections 2.2.1, 3.2 and 4.1, etc.) were the section labels in the manuscript.

Satellite-Based Tmin Estimation
A Tmin estimation model had been established based on MYD11A1 nighttime LST averages for 3 × 3 pixel windows and DEM data, employing the standard lapse rate of air temperature with elevation [32]. Due to limitations encountered from cloud contamination or other reasons, daily inconsistent LST pixels labeled by QC values of 2 or 3 were gapfilled. In this study, a clear sky MODIS LST reconstruction method (RSDAST) [34] was used. RSDAST requires only one parameter (the MODIS LST product itself) to reconstruct missing LST values based on the assumption that differences in LST between nearby pixels are relatively stable during a short time period [34]. Validation on remotely sensed Tmin estimation with 145,210 ground-based Tmin observations showed that the accuracy of this method was acceptable with a relatively high coefficient of determination (R 2 = 0.841), low root mean square error (RMSE = 2.15 • C) and mean absolute error (MAE = 1.66 • C), and reasonable normalized RMSE (NRMSE = 25.4%) and Nash-Sutcliffe model efficiency (EF = 0.12) [32].
Spatially continuous Tmin patterns at a daily scale were produced based on remotely sensed MODIS LST and ASTER GDEM products by employing the previously established Tmin estimation model and RSDAST method. Shrub classification in PFT scheme was extracted as tea planting area with MODIS land cover type products.

Spatiotemporal Characteristics of SFD
According to the Chinese Meteorological Industry Standard, daily Tmin of 4 • C was selected as the trigger temperature of SFD for tea plants in the MLRYR [17]. In this paper, the temporal dynamics and spatial characteristics of spring frost damage to tea plants were analyzed in the MLRYR based on the estimated Tmin below 4 • C.
For the temporal dynamics of SFD, the analysis focused on two temporal scales: (i) Annual dynamics. The SFD days for each pixel in the tea planting areas over the MLRYR were counted from March to April in 2003-2020, and then the average days of SFD in the total tea planting pixels were calculated for each year. (ii) Daily dynamics. The SFD pixels for each day in the tea planting areas over the MLRYR were counted from March to April in 2003-2020, and daily pixels of SFD were averaged from 2003 to 2020. Meanwhile, the SFD days for each station and the SFD stations for each day at 163 groundbased meteorological stations over the MLRYR were also counted from March to April in 2003-2020, in order to further illustrate the consistency of Tmin between ground-based meteorological observations and remotely sensed estimations at the two temporal scales.
For the spatial characteristics of SFD, the SFD days for each pixel in the tea planting areas over the MLRYR were counted from March to April in 2003-2020, and then the average days of SFD were calculated in 2003-2020. Additionally, annual SFD days at 163 ground-based meteorological stations were also used to analyze the spatial characteristics of SFD during the past two decades at the multi-year average scale and the individual year scale.

Multidimensional Risk of SFD
Three scales of risks of SFD, including temporal, spatial (latitude and longitude), and terrain scales, were analyzed. For temporal risk, daily pixel-based SFD percentage from March to April was calculated as the ratio of daily SFD pixels in every year to the total tea planting pixels. For spatial and terrain risks, the pixel-based SFD frequencies were averaged by grouping different latitude, longitude, and elevation over the total tea planting pixels in the MLRYR from March to April in 2003-2020.

Mann-Kendall Test
The Mann-Kendall (M-K) statistical test [35,36] is a rank-based nonparametric test, and it has been widely used to detect monotonic trends in series of environmental, climate, or hydrological data [37][38][39].
The sequential version of the M-K test [40] was used to test the start time of a trend within the sample X 1 , X 2 , . . . . . . , X n from a set of random variable X based on rank series of progressive and retrograde rows of the sample. The sequential values of the test statistic U(t) were calculated with the mean (E(t)) and variance (Var t j ) of the test statistic (t j ) as done by Gerstengarbe and Werner [41] and Zhao et al. [42]: where t j is the cumulative number in which x i is greater than x j (1 ≤ j ≤ i), U(t) is the forward sequence, and is denoted as UF(t). The backward sequence U (t) was calculated using the same equation but in the reverse data series, and is denoted as UB(t) in this paper. When the test statistic UF is greater than 0, the monotonic trend is increasing, while a negative UF value indicates a decreasing trend [39]. Additionally, the M-K test can be used to detect whether the trend is statistically significant or not by checking the critical value of UF from the standard normal table. For instance, the critical limit is ±1.96 when the significance level is 0.05 [38], and the limit is ±2.56 when the significance level is 0.01. The M-K test can also be used to detect cross points (also called mutation points) in temporal trends when two statistics, UF and UB, intersect within the range of ±1.96 [43,44]. The inter-annual SFD days for tea plants from 2003 to 2020 were analyzed using regional average SFD days over total tea planting pixels in the MLRYR (Figure 3). During the past two decades, SFD days exhibited a fluctuating but significantly (p < 0.05) decreasing trend (2.7 days per ten years). The maximum SFD days appeared in 2010 with a value of 16.5 d, and the minimum SFD days occurred in 2013 with a value of 7.5 d. Furthermore, annual standard errors (SE) and 95% confidence intervals (CI) of SFD days were also calculated over total tea planting pixels. The results showed that annual SEs were between 0.03 and 0.04 d, and 95% CIs were between 7.40 and 16.56 d (not show on Figure 3). The inter-annual SFD days for tea plants from 2003 to 2020 obtained from regional pixels were plotted against ground-based average SFD days from 163 meteorological stations in the MLRYR (Figure 4). From the scatter plot, we can see that annual SFD days averaged from all regional tea-planting pixels were systematically consistent with the annual averages obtained from ground-based meteorological observations during the 18-year study period. The linear regression coefficient and the coefficient of determination (R 2 ) were 0.95 and 0.58, respectively, and were significant at p < 0.01. However, the linear regression intercept was greater than 0 (5.4), indicating that systematic deviation still existed.

Daily Dynamics
Daily dynamics of regional pixels and meteorological stations affected by SFD from March to April are shown in Figure 5. The regional pixels and meteorological stations affected by SFD, to some extent, can represent SFD coverage. Generally, SFD coverage decreased with increasing date at a rate of 244.56 pixels/day, which was equivalent to 2.45 × 10 4 ha/day. The highest coverage occurred in the first half of March, followed by a sharp decrease in coverage during the next 30 days, and then the coverage dropped to the lowest values at the end of April. The largest SFD coverage denoted with regional pixels appeared on March 6 with 17,862 pixels, accounting for 43.4% of the entire tea planting area in the MLRYR. Correspondingly, the maximum number of SFD stations occurred on March 4, and this date was very close to the date of the largest number of regional pixels. The number of SFD stations on March 4 was 69, accounting for 42.3% of the total 163 groundbased meteorological stations. The two percentages (regional pixels and meteorological stations) were similar, indicating that the area influenced by SFD was basically identical by the two assessments. The consistency between the two datasets of regional pixels and meteorological stations was further checked by linear regression. The regression results revealed that the two datasets were highly correlated, with R 2 of 0.81 for the regression (significant at p < 0.01) (in the small inset panel in Figure 5).
The maximum, minimum, and average SFD pixels during 2003-2020 were also shown in Figure 5. Generally, all three SFD pixels fluctuated over time, but declined with increasing date. Almost all of the daily maximum SFD pixels were greater than 20,000 before March 25, indicating that there was at least one year in which more than half of the tea planting areas (the number of pixels in the tea planting areas is 41,175) in the MLRYR suffered from SFD threats. After that date, although the maximum SFD pixels was gradually decreasing, the SFD coverage was still increasing from time to time. For example, the maximum SFD pixels dropped to 5000 around April 20. However, the SFD pixels increased again to more than 12,000 during the last several days of April, and even reached near 28,000 on April 27. The minimum SFD pixels was relatively low with values of less than 3300 from March 1 to April 30, indicating that the SFD in some years was almost none and the tea benefits may have been great.   Overall, the impact of spring frost was scattered across almost all shrub regions in the MLRYR. The regions with the greatest SFD days, denoted with the dark blue color, were distributed mainly in the northwestern part of MLRYR. The average SFD days during the past 20 years in northwestern MLRYR were not less than 16 d, accounting for more than 26% of the total days from March to April. There were many days in this area where there was severe SFD. Areas with few SFD days, expressed with red color for 5 d or less, were primarily located in the southern region of the MLRYR. The areas with SFD days between 6 and 10 d were widely distributed, and included parts of north-central Hunan and Jiangxi Provinces, parts of southern Hubei and Anhui Provinces, and the border areas of Jiangsu and Zhejiang Provinces. Those regions are also the major tea planting areas in the MLRYR. where there were red and purple pixels. Additionally, the smaller circles were also located in the northwestern part of MLRYR, where the SFD days were 5 d less than the satellitebased results. The possible reason for this slight difference between satellite-based and ground-based SFD days in northwestern MLRYR was that the altitude at the meteorological station was relatively low compared with the tea gardens, and the SFD days determined from satellite-based Tmin estimations at the tea garden site altitude and the SFD days determined from ground-based Tmin observations at the meteorological observation site altitude were naturally different.

Annual SFD Days
Spatial characteristics of annual SFD days, denoted by both continuous pixel-based Tmin in different colors and by discrete ground-based Tmin in different size open circles, in 2003-2020 were individually displayed for each year in Figure 8. As stated earlier, the northeastern region with larger open circles at meteorological stations had almost no tea planting, so these areas can be ignored. For other regions with many tea gardens, there were generally more SFD days in the northern region, and fewer days in the southern MLRYR. The areas with larger circles basically coincided with the scattered pixels in blue colors. On the other hand, the regions with smaller circles coincided with the red and purple pixels. with ground-based minimum temperature (Tmin) observations, respectively. Blue frames around panels indicate years with large areas covered by more SFD days, red frames around panels indicate years with large areas covered by fewer SFD days, and no frames indicate years with nearly even areas covered by more SFD days and fewer SFD days.
The number of years with large, medium, and small areas covered by more SFD days were relatively equal in 2003-2020 (five, six, and seven years, respectively). On one hand, the years covered mostly by the blue and dark blue pixels and the large open-circle areas included 2003, 2010, 2011, 2016, and 2017 (panels framed with blue squares in Figure 8). The areas with blue and dark blue pixels signifying more SFD days and with large circles were mainly located in the northern and central MLRYR, where annual SFD days were 11 days or greater. The regions with five or fewer SFD days (red pixels and small circles) were located in the southern parts of Hunan and Jiangxi Provinces. On the other hand, the number of years having larger red-colored areas was seven, including 2006,2008,2013,2014, and 2018-2020 (panels framed with red squares in Figure 8). The red and purple colored areas almost covered the entire Jiangxi and Anhui Provinces, most of Hunan, Jiangsu, and Zhejiang Provinces, and the central and eastern parts of Hubei Province. The blue and dark blue colored areas were only poorly concentrated in the northwestern parts of the MLRYR, and occasionally scattered in the other regions not mentioned above. Finally, the number of years in which the areas covered by more and fewer SFD days were nearly equivalent was six, including 2004,2005,2007,2009,2012, and 2015 (without frames in Figure 8). The areas with dark blue pixels and large circles mainly concentrated in the north of the MLRYR, the regions with blue pixels and medium circles generally distributed in the central parts of the MLRYR, and the areas with warm tones and small circles primarily clustered in the south of the MLRYR. Generally speaking, the number of years with smaller SFD areas and fewer SFD days accounted for the majority during the past twenty years, especially from 2011 to 2020.

Temporal Risk
Daily risk of SFD, denoted by the pixel percentage affected by SFD, over all tea planting areas in the MLRYR from March to April in 2003-2020 is presented in Figure 9. Among the 18 years studied, there were no years in which the average pixel-based percentage of SFD exceeded 40% (see rightmost column of Figure 9), and the number of years in which this proportion was no more than 20% or between 20-40% was nine years. Furthermore, the years with the average pixel-based percentage of no more than 20% were mainly concentrated in the past 10 years. By focusing on the daily average pixel-based percentages of SFD, we can see that the days with percentages greater than 40% were mainly concentrated in the first ten days of March, and the smaller percentages of less than 20% generally occurred in April (see topmost row of Figure 9). Daily pixel percentages of SFD from March to April also decreased with increasing date at a rate of −0.59%/d. The total of 1098 (61 days by 18 years) daily pixel-based percentages also revealed that the higher percentages were concentrated in the first twenty days of March in the earliest years, and the smaller percentages normally appeared in the last twenty days of April in the past ten years. Even so, in mid-to-late April, there were still some years when the pixel percentage of SFD exceeded 40%, and this percentage was unexpectedly up to 68% in 2007. It should be pointed out that the percentage in the first ten days of April was not too high in most years. However, the number of days with percentages exceeding 60% was one day in 2009, one day in 2010, and two days in 2020. The maximum percentage occurred on April 4, 2020 with a value of 86%, indicating that more than four fifths of the tea planting area in MLRYR was affected by SFD. Therefore, the temporal risk of SFD for tea planting areas has nearly always existed, and it cannot be ignored.

Spatial Risk
Spatial risk of pixel-based SFD, quantified by the frequencies of occurrence of pixelbased SFD at different latitudes and longitudes (0.1 • intervals) averaged over all of the tea planting areas in the MLRYR from March to April in 2003-2020 is shown in Figure 10. When latitude was less than 33.2 • N, the average frequency for SFD was raised with increasing latitude (Figure 10a). In areas with relatively higher latitudes, such as at latitudes greater than 33.2 • N, the average SFD frequencies fluctuated drastically due to the smaller pixel numbers of tea gardens.
From the perspective of longitude, the average SFD frequencies in the westernmost part of the study area (the area with lower longitude), were relatively higher compared with other longitudes. SFD frequencies gradually decreased with increasing longitude, and then SFD values were basically stable for longitude greater than about 112 • E (Figure 10b). The areas at relatively smaller longitudes had higher SFD frequencies, and this observation may be related to the long distance from the ocean or the higher altitude.

Terrain Risk
Terrain risk of pixel-based SFD, quantified by the frequencies of occurrence of pixelbased SFD at different altitudes (100 m intervals), averaged over all of the tea planting areas in the MLRYR from March to April in 2003-2020 is shown in Figure 11. The minimum values of all five statistical variances, including the maximum, the minimum, the median, the first quartile, and the third quartile, occurred for the regions where altitude was less than 100 m. The medians of all SFD frequencies for all altitudes were greater than 14.3%, and the maximum median value of 68.4% occurred at 1600 m. The increase in SFD frequency with increasing altitude was statistically significant. Only at the two highest altitudes did the SFD frequencies decrease irregularly. Furthermore, both the difference between the maximum and the minimum frequencies and the difference between the third quartile and the first quartile frequencies were relatively stable with altitude, and the differences were around 18% and 6.6%, respectively, except for the highest three altitudes. The narrow ranges within the same altitude interval indicated that the inter-annual fluctuations among different altitudes were very small. Notably, the change in SFD frequency (box plots, left y-axis) and SFD pixel number (black line with open circles, right y-axis) with altitude presented opposite patterns. Areas with higher altitudes had a higher SFD frequency, but fewer tea plantations were located in these areas. In contrast, in areas where there were more tea plantations, SFD frequency was relatively low.

Advantage for Satellite-Based Tmin Estimation
Daily minimum temperature is the critical trigger threshold of SFD for crops (winter wheat [45], winter sugar beet [46], potato [47]); orchards (apple and cherry [11], tea [14]); and forests (beech [10,48,49], spruce [11]). The spatiotemporal distribution of Tmin, therefore, needs to be accurately identified in order to know whether SFD has occurred and what the extent of the impacts are. However, ground-based Tmin observations are discontinuous and heterogeneous because of limitations in the weather station network layout. The discrete distribution of weather stations will result in underestimation or overestimation for spatial SFD risk, especially in areas where the meteorological stations are themselves very sparse [21], and in areas where elevation changes drastically. Due to the availability of LST data in recent years for regions with inadequate or inaccessible weather station networks, the retrieval of minimum air temperature from remotely sensed LST has been a popular method for estimating Tmin, and good relationships have been reported between nighttime LST and Tmin [21][22][23]. By employing the standard lapse rate of air temperature with elevation, Wang et al. established a Tmin estimation model based on MODIS nighttime reconstructed LST values and DEM data [32]. This simple estimation model effectively solved the problem of spatially discontinuous and heterogeneous weather station locations, and was further confirmed to be a useful method for monitoring and evaluating SFD risk for tea plants in the MLRYR of China.

Rationality of Risk Distribution
SFD is a well-known worldwide risk for overwintering crops, tree species, and some early spring orchards in temperate climate regions [8,50,51] and Mediterranean areas [49,52] at northern latitudes, and also for tea plantations in humid areas around the equator [14]. The distribution characteristics of SFD risk were systematically analyzed at three scales, including temporal scale, spatial scale, and terrain scale, over the MLRYR in China. The results of our study demonstrated that the period with higher SFD risk for tea plants occurred primarily in early spring over the past two decades. Even so, the SFD risk during late spring (e.g., April) deserves to be given more attention due to occasionally drastic cooling events. Recent research in the Southern Apennine Mountains found that lower temperatures in May in 2019 lengthened the duration of beech leaf development (from dormant bud to leaves completely unfolded) by nine days, and the date of unfolded leaves was even delayed 11 days later compared with the number of days in 2018 [52].
Spatially, the regions with higher SFD risk in the MLRYR were concentrated in the northwestern parts of the study area where altitude was also relatively higher. The results found in this study have also been observed in other studies. For example, Laughlin found that minimum temperatures were lower at higher elevations in the absence of cloud cover and disturbances by wind speed [53]. Kotikot et al. also found that higher elevations increased the odds of a frost occurrence, even though there was only a 400 m elevation range [54]. Additionally, Laughlin and Kalma examined watershed area as a factor affecting frost occurrence, and found that it could explain 65% of the unexplained variation in minimum temperatures [55]. This result may explain the higher SFD risk appearing in the westernmost part of the study area. Kotikot et al. proved that slope aspect was an important factor for frost occurrence [54]. However, different researchers have come to inconsistent conclusions over different regions. Kotikot et al. thought that the odds of frost occurrence on a western slope in equatorial regions were reduced by a factor of 0.12 compared with an eastern slope [54]. Longstroth reported that northern facing slopes were more susceptible to frosts than southern facing slopes in North America [56]. Gurskaya et al. found that northeastern facing slopes were the most frost-prone in the southern Ural Mountains [57]. Therefore, more work should be done in the future to further evaluate SFD risk by considering more potential factors, such as land cover types, distance to water exposure, slope, aspect, etc.

Impacts of Spring Frost Types
Spring frost is a catastrophic weather event caused by the intrusion of cold air. Frosts can be classified as advection frost, radiation frost, or mixed frost [58]. Advection frost occurs when cold air accumulated at high latitudes moves southward, and consequently drops temperature drastically and is followed by strong winds. Radiation frost is caused by strong radiational cooling on the ground surface under the control of cold high pressure that produces clear and windless conditions. Radiation frosts are often very localized and site-specific [59]. They typically occur during springtime when there is a combination of cloudless nights with little or no wind, low humidity, and low temperature. Mixed frost is a combination frost arising from cold air advection from strong winds and radiational cooling. Different frost types may have different impacts on vegetation, depending on where the vegetation grows. Advection frost is mainly found on vertical exposed surfaces, and radiation frost usually occurs in valleys since cold air is heavier than warm air [60]. No matter the type of spring frost, once it occurs, it will damage newly formed buds and flowers, destroy young fruit and crops, and even devastate full trees, thereby bringing huge economic losses to local farmers. It is essential to have prior knowledge of species-specific responses to frost events [60]. The good news is that there are ways to control spring frosts and to mitigate the subsequent damage to plants. Common methods used to protect plantations from frost include smokers, wind machines, and propane burners. Naturally, there are pros and cons to each method. The proper choice depends on many factors such as topography, cost of equipment, and labor requirements (Save money on frost protection with sprinkler irrigation I Netafim (netafimuk.com. Last visited on 20 May 2021)). Effective adaptations may benefit from genetic developments and/or phenotypic plasticity, and the latter is likely to be more important [61].

Potential Migration of Risk under Climate Warming
Climate warming has been widely accepted. According to climate model projections, the global average air temperature is likely to increase by 0.3-4.8 • C by the end of this century [2]. Temperature is commonly regarded as one of the primary drivers of spring vegetation phenology [62], especially for perennial trees at middle and high latitudes of the Northern Hemisphere [63], and various vegetation types (deciduous needle-leaved forest, alpine and sub-alpine meadow, plain grassland, desert grassland) in Xinjiang in the center of the Eurasian continent [62]. A notable longitudinal change pattern for natural ecosystems has been found, with considerable increases in temperature sensitivity from inland areas to most coastal areas over the North Hemisphere [64]. On one hand, warmer conditions in winter might reduce the effective chilling accumulation, and therefore delay spring phenology [65]. On the other hand, temperature increases in early spring could advance the fulfillment of heat requirements and bring earlier spring phenology [66]. However, climate change might further expose young buds to unstable weather environments, leading to an increased risk of spring frost damage. Ideally, if spring phenology advances more slowly than late spring frost threats, SFD risk may be reduced. Recent research demonstrated that delayed blossom of apple trees was observed from 2001 to 2016 in Shaanxi, China, and late spring frosts mostly occurred before tree flowering, causing only minor frost damage to apple trees [65].

Uncertainties and Limitations
The MLRYR region, including six provinces and one municipality, was analyzed as a single area. The total land area in the MLRYR exceeds 1,010,000 km 2 , extending from 24 • to 35 • N latitude and from 109 • to 122 • E longitude. The trigger threshold of spring frost for tea plants was set as 4 • C based on an experimental conclusion using variety Longjing 43 as an example [67]. In fact, there are many varieties of tea plants varying over a large spatial extent, and their tolerances for frost damage, consequently, may be different. Therefore, a detailed study of trigger thresholds for frost damage should be done to precisely evaluate the SFD risk for different varieties of tea plants. Moreover, the dates of tea budding and leaf unfolding may come earlier under climate warming, and a phenology-based approach for estimating the tea picking period should be studied in order to accurately assess crop-specific frost risk.

Conclusions
Results in this study revealed that both SFD days at the annual scale and SFD areas at the daily scale exhibited a decreasing trend of 2.7 days/decade and 2.45 × 10 4 ha/day, respectively. Spatially, the regions with the greatest SFD days were mainly distributed in the northwestern part of MLRYR, and the areas with few SFD days were primarily located in the south of MLRYR. Additionally, SFD risks at three scales, including temporal, spatial, and terrain scales, were also evaluated. Generally, higher SFD risks were found in the first twenty days of March at the temporal scale, in the northwestern region at the spatial scale, and in the relatively higher elevation areas at the terrain scale.
This study further confirmed that estimated Tmin based on MYD11A1 nighttime products and DEM is a useful way to map and evaluate threats and potential risks of spring frost damage to tea plantations in the MLRYR China. Our findings have important implications for government agencies to reasonably determine layouts for new tea gardens at higher elevations, for agricultural insurance agencies to draft effective tea insurance strategies at different regions, and for tea farmers to give greater attention to potential frost damage and to take appropriate frost prevention measures to reduce losses.  Data Availability Statement: All data used during the study area available at http://data.cma.cn/ site/index.html, https://earthdata.nasa.gov/, and https://doi.org/10.5067/ASTER/ASTGTM_NC. 003 (accessed on 25 April 2021).

Conflicts of Interest:
The authors declare no conflict of interest.