Spatial Heterogeneity in the Response of Winter Wheat Yield to Meteorological Dryness/Wetness Variations in Henan Province, China

: Knowledge of the responses of winter wheat yield to meteorological dryness/wetness variations is crucial for reducing yield losses in Henan province, China’s largest winter wheat production region, under the background of climate change. Data on climate, yield and atmospheric circulation indices were collected from 1987 to 2017, and monthly self-calibrating Palmer drought severity index (sc-PDSI) values were calculated during the winter wheat growing season. The main results were as follows: (1) Henan could be partitioned into four sub-regions, namely, western, central-western, central-northern and eastern regions, based on the evolution characteristics of the time series of winter wheat yield in 17 cities during the period of 1988–2017. Among them, winter wheat yield was high and stable in the central-northern and eastern regions, with a remarkable increasing trend ( p < 0.05). (2) The sc-PDSI in February had significantly positive impacts on climate-driven winter wheat yield in the western and central-western regions ( p < 0.05), while the sc-PDSI in December and the sc-PDSI in May had significantly negative impacts on climate-driven winter wheat yield in the central-northern and eastern regions, respectively ( p < 0.05). (3) There were time-lag relationships between the sc-PDSI for a specific month and the atmospheric circulation indices in the four sub-regions. Furthermore, we constructed multifactorial models based on selected atmospheric circulation indices, and they had the ability to simulate the sc-PDSI for a specific month in the four sub-regions. These findings will provide scientific references for meteorological dryness/wetness monitoring and risk assessments of winter wheat production.


Introduction
Global warming has led to an ever-intensifying water cycle, and in particular, increasing rates of evaporation and precipitation, which have also changed the occurrence and severity of meteorological dryness/wetness events in different regions of the world [1][2][3][4].Crop farming has been the backbone of human society, providing food and other basic needs for people [5].With the rapid development of urbanization and industrialization, a large amount of cultivated land has been occupied, so improving crop yield per unit area has proven to be an important approach to boost food production [6,7].However, the weight of observational evidence indicates that unfavorable meteorological dryness/wetness events Agronomy 2024, 14, 817 2 of 15 during the crop growing season can lead to reduced crop yields and even threaten food security [8][9][10].Therefore, analyzing the spatial heterogeneity of yield variability and its response to meteorological dryness/wetness conditions is of great importance in evaluating regional yield vulnerability and proposing adaptation strategies under global warming.
With the increasing importance of long time-series data, time series clustering methods have become widespread in the field of hydrometeorology [11,12]; they divide a large amount of data into smaller groups according to their similarity so that the data belonging to the same group are as similar as possible.Traditional clustering methods designed for static data are susceptible to noise when dealing with the time series clustering problem due to the high dimensionality of the time series [13].Under these circumstances, the characteristic-based time series clustering method was proposed [14], which clusters based on structural characteristics extracted from the time series rather than using a distance measure to cluster point values.The characteristic measures used to reflect the time series are derived from statistical operations, such as statistical, time-domain and frequency-domain characteristics, and they are then input into arbitrary clustering algorithms.For example, Garcia-Gutierrez et al. (2022) [15] used a K-means clustering algorithm to regionalize solar radiation in Spain based on statistical parameters computed from clear-sky index series.As a typical example of time series data, the regionalization of crop yield for a given region not only analyzes spatial variations in yield but also provides fine-grained strategies to cope with unfavorable meteorological dryness/wetness conditions [16].Despite the widespread interest in the mean values of and trends in crop yield [17], the statistical, time-domain and frequency-domain characteristics of yield data have rarely been systematically examined in previous studies.
To understand the response of crop yield to meteorological dryness/wetness variations, several meteorological dryness/wetness indices have been proposed in recent decades [3,18,19], such as the standardized precipitation index (SPI), standardized precipitation evapotranspiration index (SPEI) and self-calibrating Palmer drought severity index (sc-PDSI).Both SPI and SPEI have similar multi-timescale characteristics, with the former considering only precipitation and the latter considering precipitation and potential evapotranspiration [20].In contrast to the SPI and SPEI, the sc-PDSI is more physical, as it takes into account the water balance using a two-layer bucket model and uses self-calibrating processes that can adjust the original PDSI to the local climate [21].Recent studies have suggested that the sc-PDSI demonstrates good performance in monitoring meteorological dryness/wetness variations in China, and studies have also been carried out to show that climate-yield relationships vary between regions and crop types [22,23].Taking China as an example, the sc-PDSI showed positive relationships with maize yield anomalies in Northeast China and wheat in North, South and Northwest China, but some regions (rice in East and Southeast China) had negative correlations [23].These existing studies have focused mainly on relationships by considering the entire crop-growing season as a whole [24,25].Few studies have examined seasonal differences in yield responses to meteorological dryness/wetness conditions, particularly the key month affecting yield and the corresponding key meteorological drought/wetness index (KMDWI).Considering that atmospheric circulations can be regarded as important premonitory influencing signals in the analysis of meteorological dryness/wetness variations [26,27], it is necessary to investigate the relationships between the KMDWI and atmospheric circulation indices in order to construct an appropriate simulation model, which may be a valuable approach to sustainable crop production and agricultural water management.
Wheat is an ancient and important grain crop.China is the world's largest winter wheat producer, accounting for nearly 20% of the global wheat yield [28].As the largest winter wheat production region in China, Henan province accounts for nearly 25% of the winter wheat yield in China [29].However, an uneven spatio-temporal distribution of precipitation has severely affected sustainable winter wheat production in Henan [30].Although drought and flood events frequently occur in Henan [31], the spatial heterogeneity of yield variability and its response to meteorological dryness/wetness conditions are still unclear in this region.Therefore, the main research objectives are (1) to understand the spatial difference in winter wheat yield using a characteristic-based time series clustering analysis and (2) to identify the KMDWIs that affect yield fluctuations in different subregions and evaluate the potential impacts of atmospheric circulation indices.The detailed results are of great significance for meteorological dryness/wetness monitoring and risk assessments of winter wheat production.

Study Region
Located in central China (Figure 1), Henan province covers an area of 167,000 km 2 .It consists of 18 cities, i.e., Zhengzhou, Kaifeng, Luoyang, Pingdingshan, Anyang, Hebi, Xinxiang, Jiaozuo, Puyang, Xuchang, Luohe, Sanmenxia, Shangqiu, Zhoukou, Zhumadian, Nanyang, Xinyang and Jiyuan, and it has a wide altitude ranging from 16 to 2272 m.Henan is located mainly in the warm temperate zone, with four distinct seasons and frequent meteorological disasters [31].The annual average temperature ranges from 10.5 • C in the north to 16.7 • C in the south, and the annual precipitation ranges from 407.7 mm to 1295.8 mm, mainly concentrated between June and August [16].As a major grainproducing province in China, the output of major agricultural products such as grain, cotton and oil in Henan ranks at the forefront of the country.Among them, winter wheat is the first grain crop in Henan, with 1/16 of the nation's arable land producing 1/4 of the winter wheat in China [29].
Although drought and flood events frequently occur in Henan [31], the spatial hetero neity of yield variability and its response to meteorological dryness/wetness conditi are still unclear in this region.Therefore, the main research objectives are (1) to underst the spatial difference in winter wheat yield using a characteristic-based time series c tering analysis and (2) to identify the KMDWIs that affect yield fluctuations in differ sub-regions and evaluate the potential impacts of atmospheric circulation indices.The tailed results are of great significance for meteorological dryness/wetness monitoring risk assessments of winter wheat production.

Study Region
Located in central China (Figure 1), Henan province covers an area of 167,000 km consists of 18 cities, i.e., Zhengzhou, Kaifeng, Luoyang, Pingdingshan, Anyang, H Xinxiang, Jiaozuo, Puyang, Xuchang, Luohe, Sanmenxia, Shangqiu, Zhoukou, Z madian, Nanyang, Xinyang and Jiyuan, and it has a wide altitude ranging from 16 to 2 m.Henan is located mainly in the warm temperate zone, with four distinct seasons frequent meteorological disasters [31].The annual average temperature ranges from 1 °C in the north to 16.7 °C in the south, and the annual precipitation ranges from 407.7 m to 1295.8 mm, mainly concentrated between June and August [16].As a major grain-p ducing province in China, the output of major agricultural products such as grain, cott and oil in Henan ranks at the forefront of the country.Among them, winter wheat is first grain crop in Henan, with 1/16 of the nation's arable land producing 1/4 of the win wheat in China [29].

Datasets
The homogenized datasets of monthly temperature, precipitation, atmospheric p sure, relative humidity, sunshine and wind speed from 116 stations in Henan during period of 1987-2017 were obtained from the National Climate Center of China.Data

Datasets
The homogenized datasets of monthly temperature, precipitation, atmospheric pressure, relative humidity, sunshine and wind speed from 116 stations in Henan during the period of 1987-2017 were obtained from the National Climate Center of China.Data on annual winter wheat yield (AWWY, unit: kg/hm 2 ) in 17 cities in Henan, with the exception of Jiyuan, during the period of 1987-2017 were obtained from the Henan Province Bureau of Statistics.Monthly atmospheric circulation indices that may drive the regional climate in Henan were used to analyze the premonitory influencing signals of meteorological dryness/wetness variations, which included Niño 1+2, Niño 3, Niño 3.4, Niño 4, North Atlantic Oscillation (NAO), Arctic Oscillation (AO) and Pacific Decadal Oscillation (PDO), and they were obtained from http://www.cpc.ncep.noaa.gov/products/precip/CWlink/(accessed on 5 November 2023).All data used in this study were obtained from official sources and were subjected to rigorous quality controls.

Characteristic-Based Time Series Clustering Analysis
In this study, AWWY was selected for a given indicator, and a series of characteristicbased indices were constructed and are presented in Table 1.As shown in Table 1, we calculated the statistical characteristics of the AWWY time series, such as the minimum value (MIN), maximum value (MAX), mean value (AVE), range value (RNG), standard deviation (STD) and coefficient of variation (CV).Furthermore, a trend analysis is one of the most important measurements in studying time series data, and it can reflect their time-domain characteristics.In this study, the Z statistic (Z MK ), Sen's slope (MAG) and the Hurst index (H) were used to describe the trends and persistent patterns of the AWWY time series using the Mann-Kendall test, Sen's slope estimator and the Hurst method, respectively [32].In addition to a trend analysis, multi-time scale characteristics are also important in studying time series data.For example, the AWWY time series can be decomposed into three periodic oscillation intrinsic mode functions (IMFs) and a trend component, and this was achieved in this study using ensemble empirical mode decomposition (EEMD) [33].The variance contribution of each IMF (i.e., CIMF1, CIMF2 and CIMF3) and the trend component (CTrend) to AWWY was selected to reflect the frequency-domain characteristics of the AWWY time series.More details about the Mann-Kendall test, Sen's slope estimator, the Hurst method and EEMD are shown in the Supplementary Materials [34][35][36][37][38][39][40][41].Taking into account the potential collinearity and redundancy among these indices, a principal component analysis (PCA) and a K-means clustering analysis were used in our study to partition the 17 cities into k clusters, and the clustering result was evaluated using silhouette coefficients (SCs) [11,42].The highest SC value determined the optimum number of clusters.More details about the PCA and K-means clustering analysis are shown in the Supplementary Materials [43][44][45][46][47][48][49][50].

Calculation of sc-PDSI
The sc-PDSI is a typical meteorological dryness/wetness index based on a physical water balance model, calculated from precipitation and potential evapotranspiration, together with parameters related to the characteristics of the soil/surface at each location [51].In the calculation of the sc-PDSI, potential evapotranspiration is based on the Penman-Monteith equation, which takes into account energy balance and aerodynamic theory, rather than the Thornthwaite method [52].The sc-PDSI values range from −4 to +4, with a negative (positive) value indicating a dry (wet) period.In this study, we calculated the monthly sc-PDSI during the winter wheat growing season (from October of the previous year to May of the current year) in Henan using climate data on temperature, precipitation, relative humidity, sunshine and wind speed to reflect meteorological dryness/wetness conditions, and this was achieved using an open-access R package (https://github.com/Sibada/scPDSI,accessed on 5 November 2023).

Calculations of Climate-Driven Winter Wheat Yield (CDWWY)
Crop yield is affected by a variety of factors such as crop varieties, agricultural technology and meteorological conditions, which is a combination of trend yield and climatedriven yield [53]. Extracting the climate-driven yield is a crucial step in understanding the yield impacts of meteorological dryness/wetness variations [9].In this study, four statistical methods, namely, EEMD, linear regression (LR), logistic function (LF) and the Hodrick-Prescott (HP) filter, were used to extract CDWWY [54].
In this study, we mainly focused on the relationships between time series variables (i.e., CDWWY and sc-PDSI) during the period of 1987-2017.Therefore, we used a Pearson correlation analysis to examine the relationships between the CDWWY time series and sc-PDSI time series during the winter wheat growing season, which allowed us to select the key month that affected the yield and the corresponding KMDWI in Henan.Furthermore, a regression equation was constructed between the KMDWI and CDWWY using the following formula [42]: where std (CDWWY) and std (KMDWI) are the standardized values of the KMDWI and CDWWY, respectively.Coefficient a is the sensitivity of the yield to the KMDWI, and it is multiplied by the trend of the KMDWI for a given period to determine the impact of the KMDWI trend on yield [42].

Relationships between the KMDWI and Atmospheric Circulation Indices
Atmospheric circulation indices are considered important premonitory influencing signals in the analysis of meteorological dryness/wetness variations [26,27].A time-lag correlation analysis was performed to examine the relationships between the KMDWI and seven atmospheric circulation indices [55,56], and a time lag of 0-11 months was selected in this study.According to the results of the correlation analysis, we then selected the atmospheric circulation indices with the highest correlation with the corresponding KMDWI for each month as the premonitory influencing signals.Finally, we used the multiple linear regression method to construct an empirical KMDWI simulation model based on the atmospheric circulation indices selected in this study.Two indicators were used to evaluate the performance of the multiple linear regression model: the determination coefficients (R 2 ) and root-mean-square error (RMSE) [57].More details about the time-lag correlation analysis and multiple linear regression method are shown in the Supplementary Materials.

AWWY-Based Regionalization in Henan
As described in Section 2.3.1, the values of 13 indices in 17 cities in Henan were first calculated.Then, the principal components of 13 columns (13 indices) × 17 rows (17 cities) were analyzed using the PCA method (Table 2).To reduce the dimensionality of the data, we selected the first three principal components (PCs) because their eigenvalues were greater than 1, and their cumulative variances exceeded 89%.As shown in Table 2, PC1 had higher loadings on five indices, namely, Z MK , H, CIMF1, CIMF3 and Ctrend, representing the persistent trend of the AWWY time series.PC2 had higher loadings on four indices, namely, RNG, STD, CV and MAG, representing the magnitude of the change in AWWY.PC3 had higher loadings on three indices, namely, MIN, MAX and AVE, representing the high (or low) levels of AWWY.According to the first three PCs, we obtained their corresponding score sequences.A matrix of 3 columns (3 score sequences) × 17 rows (17 cities) was clustered using the K-means method.As shown in Figure 2, the SC values varied for different numbers of clusters.The highest SC value occurred when the number of clusters reached four.Therefore, 17 cities in Henan could be partitioned into four sub-regions characterized by different AWWY variations, i.e., Sanmenxia and Luoyang in the western region (Region I); Zhengzhou, Pingdingshan and Nanyang in the central-western region (Region II); Xuchang, Kaifeng, Jiaozuo, Xinxiang, Hebi, Puyang and Anyang in the central-northern region (Region III); and Luohe, Shangqiu, Zhoukou, Zhumadian and Xinyang in the eastern region (Region IV) (Table S1).

Spatial Differences in AWWY in Henan
The statistical values of 13 indices characterizing the AWWY time series in the four sub-regions are shown in Table 3.As shown in Table 3, the AVE of AWWY was the highest in Region III (5770.9kg/hm 2 ), followed by Region IV (5091.3kg/hm 2 ) and it was the lowest in Region I (3536.3kg/hm 2 ).Compared with the AVE, the RNG of AWWY was the highest in Region IV (4736 kg/hm 2 ), followed by Region III (3561.2kg/hm 2 ) and it was the lowest in Region II (2543.8kg/hm 2 ).Furthermore, the Z MK , MAG and CTrend were higher in Regions III and IV than in Regions II and I.These results indicate that in Henan, Regions III and IV had a remarkable increasing trend in AWWY, as well as a high and stable yield.

Spatial Differences in AWWY in Henan
The statistical values of 13 indices characterizing the AWWY time series in the four sub-regions are shown in Table 3.As shown in Table 3, the AVE of AWWY was the highest in Region III (5770.9kg/hm 2 ), followed by Region IV (5091.3kg/hm 2 ) and it was the lowest in Region I (3536.3kg/hm 2 ).Compared with the AVE, the RNG of AWWY was the highest in Region IV (4736 kg/hm 2 ), followed by Region III (3561.2kg/hm 2 ) and it was the lowest in Region II (2543.8kg/hm 2 ).Furthermore, the ZMK, MAG and CTrend were higher in Regions III and IV than in Regions II and I.These results indicate that in Henan, Regions III and IV had a remarkable increasing trend in AWWY, as well as a high and stable yield.

Relationships between CDWWY and sc-PDSIs during the Wheat Growing Season in Henan
In this study, we employed four commonly used methods to extract CDWWY (Figure 3).Compared with the LR, LF and HP filter, the CDWWY calculated using EEMD showed a higher correlation with the sc-PDSIs during the winter wheat growing season for the four sub-regions in Henan.Taking the CDWWY calculated using EEMD as an example, it was found that the sc-PDSIs from January to February of the current year had a significant positive impact on CDWWY in Region I (p < 0.05), and the sc-PDSI in February had the largest impact on CDWWY, which was considered the KMDWI in this sub-region.Similarly, the sc-PDSI in February also had the largest impact on CDWWY in Region II (p < 0.05), and a wet February was more favorable for winter wheat production in this sub-region.In contrast to Regions I and II, the sc-PDSI in December of the previous year and the sc-PDSI in May of the current year had significant negative impacts on CDWWY in Regions III and IV, respectively (p < 0.05), which were considered the KMDWI in both sub-regions, respectively.These results indicate that meteorological dryness/wetness conditions in different key months limited winter wheat production in Regions III and IV.
was found that the sc-PDSIs from January to February of the current year had a significant positive impact on CDWWY in Region I (p < 0.05), and the sc-PDSI in February had the largest impact on CDWWY, which was considered the KMDWI in this sub-region.Similarly, the sc-PDSI in February also had the largest impact on CDWWY in Region II (p < 0.05), and a wet February was more favorable for winter wheat production in this subregion.In contrast to Regions I and II, the sc-PDSI in December of the previous year and the sc-PDSI in May of the current year had significant negative impacts on CDWWY in Regions III and IV, respectively (p < 0.05), which were considered the KMDWI in both subregions, respectively.These results indicate that meteorological dryness/wetness conditions in different key months limited winter wheat production in Regions III and IV.The responses of CDWWY to KMDWI changes in different sub-regions and periods are shown in Figure 4.If the KMDWI increased by 1% during the study period, CDWWY increased by 0.66% and 0.46% in Regions I and II, respectively, and it decreased by 0.35% and 0.30% in Regions III and IV, respectively (Figure 4a).However, limited impacts of the KMDWI trends on CDWWY were found for the four sub-regions, ranging from −0.06% to 0.05%, due to a weak KMDWI trend in each sub-region from 1988 to 2017 (Figure 4b,c).Furthermore, we classified the study period into three sub-periods, i.

Relationships between the KMDWI and Atmospheric Circulation Indices in Henan
We used a Pearson correlation analysis to examine the relationship between the KMDWI of the key month and seven atmospheric circulation indices of the 12 months previous to that key month (Figure 5).As shown in Figure 5, the most significant correlation was found in Region IV, with the highest absolute value of 0.49.Compared with this, the highest correlation coefficient in Region III was less than 0.40.Furthermore, the atmospheric circulation indices had time-lag effects on the KMDWIs in the different sub-regions.For example, NAO in December, PDO in November and AO in April of the previous year had a significant impact on the KMDWI (p < 0.05) in Region I (Figure 5a), while NAO in June and December and AO in April of the previous year had a significant impact on the KMDWI (p < 0.05) in Region II (Figure 5b).As for Regions III and IV, NAO in February

Relationships between the KMDWI and Atmospheric Circulation Indices in Henan
We used a Pearson correlation analysis to examine the relationship between the KMDWI of the key month and seven atmospheric circulation indices of the 12 months previous to that key month (Figure 5).As shown in Figure 5, the most significant correlation was found in Region IV, with the highest absolute value of 0.49.Compared with this, the highest correlation coefficient in Region III was less than 0.40.Furthermore, the atmospheric circulation indices had time-lag effects on the KMDWIs in the different sub-regions.For example, NAO in December, PDO in November and AO in April of the previous year had a significant impact on the KMDWI (p < 0.05) in Region I (Figure 5a), while NAO in June and December and AO in April of the previous year had a significant impact on the KMDWI (p < 0.05) in Region II (Figure 5b).As for Regions III and IV, NAO in February and AO in August of the previous year had a significant impact on the KMDWI (p < 0.05) in Region III (Figure 5c), while the atmospheric circulation indices that affected the KMDWI in Region IV (p < 0.05) mainly included Niño 1+2, Niño 3 and Niño 3.4 from June to November of the previous year, as well as NAO in June of the previous year (Figure 5d).
Given the results of the correlation analysis described above, we selected the atmospheric circulation indices with the highest correlation with the corresponding KMDWI for each month as the premonitory influencing signals in different sub-regions.Then, we constructed regression models between the KMDWI and selected atmospheric circulation indices in each sub-region.As shown in Table 4, the multifactorial models that we constructed had the ability to simulate the KMDWIs for the four sub-regions.Their R 2 values were in the range of 0.49 to 0.65, and the corresponding RMSE values were in the range of 0.58 to 0.72 in the different sub-regions.
Agronomy 2024, 14, x FOR PEER REVIEW 10 of 16 and AO in August of the previous year had a significant impact on the KMDWI (p < 0.05) in Region III (Figure 5c), while the atmospheric circulation indices that affected the KMDWI in Region IV (p < 0.05) mainly included Niño 1+2, Niño 3 and Niño 3.4 from June to November of the previous year, as well as NAO in June of the previous year (Figure 5d).Given the results of the correlation analysis described above, we selected the atmospheric circulation indices with the highest correlation with the corresponding KMDWI for each month as the premonitory influencing signals in different sub-regions.Then, we constructed regression models between the KMDWI and selected atmospheric circulation indices in each sub-region.As shown in Table 4, the multifactorial models that we constructed had the ability to simulate the KMDWIs for the four sub-regions.Their R 2 values were in the range of 0.49 to 0.65, and the corresponding RMSE values were in the range of 0.58 to 0.72 in the different sub-regions.

Discussion
According to the AWWY time series from 17 cities, Henan province could be partitioned into four sub-regions using a characteristic-based time series clustering analysis (Figure 2).The complex geographical environment could be one of the important factors responsible for the spatial differences in AWWY.Regions I and II are surrounded by the Taihang Mountains and the Funiu Mountains (Figure 1), and the high mountain terrain restricts the use of agricultural machinery [58].Furthermore, the complex climate, poor quality of arable land and low level of irrigation in the mountainous areas have a combined effect on the production of winter wheat in both sub-regions [30,[58][59][60].However, our finding is similar to that of a risk assessment of an agrometeorological disaster in Henan, which was the result of a combination of four main factors (hazard, exposure, vulnerability and restorability) [59].Therefore, a characteristic-based time series clustering analysis could also provide a new perspective for the regionalization of crop production risk under the background of climate change due to its simple data requirements and reliable results [14].
The sc-PDSIs in February, May and December were the KMDWIs that influenced yield in the four sub-regions (Figure 3).Among them, the sc-PDSI in February had significantly positive impacts on yield in Regions I and II (p < 0.05), while the sc-PDSI in December and the sc-PDSI in May had significantly negative impacts on yield in Regions III and IV, respectively (p < 0.05).The uneven distribution of precipitation during the winter wheat growing season led to a severe water shortage in the greening stage, and the drought risk gradually decreased from western to eastern Henan [31,61].Thus, the higher the sc-PDSI in February of the current year, the higher the yields in Regions I and II.The impacts of drought risk on yield could be alleviated by maintaining soil moisture through appropriate tillage methods, improving irrigation efficiency and using drought-resistant varieties of wheat [62,63].Spring waterlogging disasters occurred mainly in the eastern parts of Henan, and a wetter climate in the maturity stage was unfavorable to spike development and grain formation [64,65], resulting in lower yields in Region IV with a higher sc-PDSI in May of the current year.Therefore, conducting farmland drainage measures in a timely manner is critical for achieving drainage and water logging, as well as coordinated regulation [66,67].Additionally, the occurrence of wheat lodging caused by heavy rainfall and high winds should also be noted [68].Chill damage in the overwintering stage was the main unfavorable factor for the growth of winter wheat in the northern parts of Henan, and more rainfall exacerbated low temperatures and low sunshine [69,70], which led to lower yields in Region III with a higher sc-PDSI in December of the previous year.Therefore, shifting the wheat planting calendar, spraying foliar fertilizers and using cold-resistant winter wheat varieties could reduce the negative effects of chill damage in the overwintering stage [71][72][73].
In addition to identifying the KMDWIs that affect yield fluctuations, we also found a time-lag relationship between the KMDWIs and atmospheric circulation indices.Although similar results have been reported in previous studies [9,55,56], our study indicates that the time-lag relationship exhibited spatial differences in the different sub-regions (Figure 5).In particular, the time lag of the same KMDWI response was different for the different atmospheric circulation indices in the four sub-regions.Furthermore, the multifactorial models used to simulate the KMDWIs for the four sub-regions were constructed by selecting suitable atmospheric circulation indices (Table 4).Compared with previous studies conducted at a larger spatial scale [74,75], the accuracy of our empirical models, together with R 2 and RMSE, could give a reliable output and be applied at the provincial scale.Furthermore, atmospheric circulation indices can be forecasted about 6 to 9 months in advance by the National Climate Center of China [76], suggesting that our empirical models obtain KMDWIs for the four sub-regions at an early stage, and they are of great significance for meteorological dryness/wetness monitoring for winter wheat production and providing some coping strategies in each sub-region [77,78].
The limitations of our study are based on two main aspects.First, we mainly focused on the response of winter wheat yield to meteorological dryness/wetness variations in Henan using statistical data from the perspective of the agricultural climate, without fully considering the effects of other yield-related factors, e.g., plant cultivation, soil category and treatments [79,80].Second, we constructed multifactorial models capable of simulating KMDWIs in four sub-regions.However, the mechanisms of the natural and anthropogenic factors, including atmospheric circulation, underlying the KMDWIs are more complex and still need to be studied in depth [81].Therefore, we will collect more data on climate, soil and crop management to understand climate-yield relationships and to improve the simulation performance of unfavorable meteorological dryness/wetness conditions by incorporating machine learning methods [82][83][84].

Conclusions
Using long-term data on climate, yield and atmospheric circulation indices, this study analyzed the response of winter wheat yield to sc-PDSI variations in Henan from 1987 to 2017.The main results were as follows: (1) A PCA and a K-means clustering analysis were used to partition the AWWY time series from 17 cities into four sub-regions.Although AWWY exhibited an increasing trend in all of the four sub-regions (p < 0.05), it was high and stable in Regions III and IV.(2) The sc-PDSI for a specific month could be considered the KMDWI affecting CDWWY in each sub-region, for example, the sc-PDSI in February of the current year for Regions I and II, the sc-PDSI in December of the previous year for Region III and the sc-PDSI in May of the current year for Region IV. (3) The atmospheric circulation indices had time-lag effects on the KMDWIs, and empirical KMDWI simulation models were constructed based on selected atmospheric circulation indices in the four sub-regions.Their R 2 values were in the range of 0.49 to 0.65, and the corresponding RMSE values were in the range of 0.58 to 0.72.

Figure 1 .
Figure 1.Geographical location of the study region and distribution of meteorological stations.

Figure 1 .
Figure 1.Geographical location of the study region and distribution of meteorological stations.

Figure 2 .
Figure 2. The clustering results of AWWY in Henan.

Figure 2 .
Figure 2. The clustering results of AWWY in Henan.

Figure 3 .Figure 3 .
Figure 3. Relationships between the climate-driven winter wheat yield (CDWWY) and the self-calibrating Palmer drought severity index (sc-PDSI) during the wheat growing season in Henan.Note: Pre-and Cur-denote the previous year and current year, respectively.A correlation coefficient value of 0.36 corresponds to the 5% significance level.
e., 1988-1999, 2000-2009 and 2010-2017, and we then found that the responses of CDWWY to KMDWI changes showed obvious interdecadal differences in Regions I and II.The sensitivities of CDWWY to KMDWI changes were the highest in 1988-1999, followed by 2000-2009 and they were the lowest in 2010-2017 in Regions I and II.Consequently, the impacts of the KMDWI trends on CDWWY gradually weakened from 1988-1999 to 2010-2017 in both sub-regions.Compared with Regions I and II, the impacts of the KMDWI trends on CDWWY were obviously weaker in Regions III and IV during the three sub-periods.changesshowed obvious interdecadal differences in Regions I and II.The sensitivities of CDWWY to KMDWI changes were the highest in 1988-1999, followed by 2000-2009 and they were the lowest in 2010-2017 in Regions I and II.Consequently, the impacts of the KMDWI trends on CDWWY gradually weakened from 1988-1999 to 2010-2017 in both sub-regions.Compared with Regions I and II, the impacts of the KMDWI trends on CDWWY were obviously weaker in Regions III and IV during the three sub-periods.

Figure 4 .
Figure 4. Responses of CDWWY to key meteorological drought/wetness index (KMDWI) changes in Henan.(a) Sensitivity of the yield to the KMDWI.(b) Trends of the KMDWI.(c) Impact of the KMDWI trend on yield.

Figure 4 .
Figure 4. Responses of CDWWY to key meteorological drought/wetness index (KMDWI) changes in Henan.(a) Sensitivity of the yield to the KMDWI.(b) Trends of the KMDWI.(c) Impact of the KMDWI trend on yield.

Figure 5 .
Figure 5. Relationships between KMDWIs and atmospheric circulation indices in Henan.Note: Preand Cur-denote the previous year and current year, respectively.A correlation coefficient value of 0.36 corresponds to the 5% significance level.

Figure 5 .
Figure 5. Relationships between KMDWIs and atmospheric circulation indices in Henan.Note: Preand Cur-denote the previous year and current year, respectively.A correlation coefficient value of 0.36 corresponds to the 5% significance level.

Table 1 .
Definitions of characteristic-based indices of annual winter wheat yield (AWWY).
Note: MIN, MAX, AVE, RNG, STD, CV, Z MK , MAG, H, CIMF1, CIMF2, CIMF3 and CTrend are abbreviations for the minimum value, maximum value, mean value, range value, standard deviation, coefficient of variation, Z statistic, Sen's slope, Hurst index and variance contribution of each IMF and a trend component, respectively.

Table 2 .
The results of the application of the principal component analysis (PCA) method to a matrix composed of 13 characteristic-based indices of AWWY in 17 cities in Henan.
Note: PC1, PC2 and PC3 are abbreviations for the first three principal components.MIN, MAX, AVE, RNG, STD, CV, Z MK , MAG, H, CIMF1, CIMF2, CIMF3 and CTrend are abbreviations for the minimum value, maximum value, mean value, range value, standard deviation, coefficient of variation, Z statistic, Sen's slope, Hurst index and variance contribution of each IMF and a trend component, respectively.

Table 3 .
Statistical values of characteristic-based indices of AWWY for four sub-regions in Henan.

Table 3 .
Statistical values of characteristic-based indices of AWWY for four sub-regions in Henan.Note: MIN, MAX, AVE, RNG, STD, CV, Z MK , MAG, H, CIMF1, CIMF2, CIMF3 and CTrend are abbreviations for the minimum value, maximum value, mean value, range value, standard deviation, coefficient of variation, Z statistic, Sen's slope, Hurst index and variance contribution of each IMF and a trend component, respectively.

Table 4 .
Regression results between KMDWIs and selected atmospheric circulation indices in Henan.

Table 4 .
Regression results between KMDWIs and selected atmospheric circulation indices in Henan.

Table 4 .
Cont.Niño 3.4 and Niño 4 reflect the Niño regions, which are areas of the tropical Pacific Ocean that are used to monitor and measure the El Niño-Southern Oscillation climate pattern, and each Niño region covers a different area of the Pacific Ocean.AO, NAO and PDO are abbreviations for North Atlantic Oscillation, Arctic Oscillation and Pacific Decadal Oscillation.