Climate-Induced Yield Losses for Winter Wheat in Henan Province, North China and Their Relationship with Circulation Anomalies

Risk analysis using climate-induced yield losses (CIYL) extracted from long-term yield data have been recognized in China, but the research focusing on the time-series characteristics of risk and the circulation signals behind yield losses still remains incomplete. To address these challenges, a case study on winter wheat production in Henan province, north China was conducted by using annual series of yield in 17 cities during 1988–2017 and monthly series of 15 types of large-scale oceanic-atmospheric circulation indices (LOACI). A comprehensive risk assessment method was established by combining the intensity, frequency, and variability of CIYL and principal component analysis (PCA). The results showed that the westernmost Henan was identified as the area of higher-risk. PCA and Mann–Kendall trend tests indicated that the southern, northern, eastern, and western areas in Henan province were classified as having different annual CIYL variations in these four sub-regions; the decreasing trend of CIYL in northern area was the most notable. Since the 2000s, a significant decline in CIYL was found in each sub-region. It should be noted that the key LOACI, which includes Tropical Northern Atlantic Index (TNA), Western Hemisphere warm pool (WHWP), and Southern oscillation index (SOI), indicated significant CIYL anomalies in some months. Furthermore, the regional yield simulation results using linear regression for the independent variables of year and various LOACI were satisfactory, with the average relative error ranging from 3.48% to 6.87%.


Introduction
As reported by Food and Agriculture Organization, global warming over the past few decades has brought frequent extreme high temperature and hydrological events, which have resulted in severe crop yield losses and reduced the income of farmers [1]. Although the improvement of crop varieties, utilization of large amount of fertilizer and effective irrigation have led to significant increases in crop production, China is particularly prone to suffer from climate disasters [2]. Therefore, the risk analysis of crop yield losses caused by climatic anomalies has been extensively conducted, mainly focusing on content, methods, indicators, and technologies [3]. Based on historical crop yield data, the extracting of yield losses from climate-induced yield (CIY) series directly reflects the comprehensive effects of unfavourable climatic conditions and even meteorological disasters; this method also effectively avoids the parametric uncertainty and data complexity present in other risk assessment models [4]. The multi-annual averaged losses rate has been especially widely applied in the risk assessment of grain crops (e.g., wheat, maize, rice) in different regions across China, owing to its simplicity and easy computation [5]. However, existing studies have mainly focused on the spatial distribution pattern of loss rates, while the exploration on time-series characteristics, including frequency, intensity, and tendency, are still limited. To date, reports identifying the inter-annual and inter-decadal variations of crop yield losses are rare.
Winter wheat is one of most important staple food crops worldwide, and China is the world's largest wheat producer and consumer [6]. Therefore, there are a large number of papers focusing on the effects of climate change on winter wheat yields in China [7]. Climatic variations are significantly impacted by large-scale oceanic and atmospheric circulation, and several studies have focused on the teleconnections between circulation patterns and wheat yields. For instance, Shuai et al. (2013) [8] pointed out that the El Niño-Southern Oscillation (ENSO) has significant impacts on wheat yields, and more rainfall in El Niño years leads to a decrease in wheat yields in the provinces of the Yangtze River basin. Xu et al. (2020) [7] indicated that the inter-decadal variations in winter wheat yields in major growing regions of China were significantly associated with the Pacific decadal oscillation (PDO), with its contribution rate to yields being approximately 11%. Chen et al. (2020) [9] found that several circulation indices, including PDO, ENSO, SOI, and NAO, had significant impacts on the inter-annual fluctuations of wheat yields in different cities of northern China. It can be concluded from previous studies that the circulation patterns should be considered as potential indication signals of wheat yield risks. Exploring the statistical relationships between climate-induced yield losses (CIYL) and different types of large-scale oceanic and atmospheric circulation indices (LOACI) may be a novel and valuable approach to yield prediction.
Henan Province is the primary producer of winter wheat in China, accounting for 25% of national production [9]. Therefore, understanding the spatial and temporal evolution of CIYL of winter wheat in Henan is important for national food security. The main objectives of this paper are: (1) to explore the spatial distribution of the risk of CIYL for wheat production in Henan by using the assessment indicators extracted from the long-term yield data in 17 cities and principal component analysis (PCA); (2) to detect regional differences in the temporal changing features of CIYL across Henan using the Mann-Kendall trend test and Pettitt's test; (3) to obtain correlations between CIYL and LOACI, calculated by considering the months during both the growing and non-growing period; and (4) to select key LOACI affecting CIYL as prediction factors for modeling actual winter yields. The detailed technology roadmap of this study is shown in Figure 1.

Climate-
climatic conditions and even meteorological disasters; this method also effectively avoids the parametric uncertainty and data complexity present in other risk assessment models [4]. The multi-annual averaged losses rate has been especially widely applied in the risk assessment of grain crops (e.g., wheat, maize, rice) in different regions across China, owing to its simplicity and easy computation [5]. However, existing studies have mainly focused on the spatial distribution pattern of loss rates, while the exploration on time-series characteristics, including frequency, intensity, and tendency, are still limited. To date, reports identifying the inter-annual and inter-decadal variations of crop yield losses are rare.
Winter wheat is one of most important staple food crops worldwide, and China is of large-scale oceanic and atmospheric circulation indices (LOACI) may be a novel and valuable approach to yield prediction.
Henan Province is the primary producer of winter wheat in China, accounting for 25% of national production [9]. Therefore, understanding the spatial and temporal evolution of CIYL of winter wheat in Henan is important for national food security. The main objectives of this paper are: (1) to explore the spatial distribution of the risk of CIYL for wheat production in Henan by using the assessment indicators extracted from the long-term yield data in 17 cities and principal component analysis (PCA); (2) to detect regional differences in the temporal changing features of CIYL across Henan using the Mann-Kendall trend test and Pettitt's test; (3) to obtain correlations between CIYL and LOACI, calculated by considering the months during both the growing and non-growing period; and (4) to select key LOACI affecting CIYL as prediction factors for modeling actual winter yields. The detailed technology roadmap of this study is shown in Figure 1.

Study Area and Data Resources
Henan province is located on the middle and lower reaches of the Yellow River in Central China (situated at 110 • 21 -116 • 39 N and 31 • 23 -36 • 22 E). As a province located in the warm temperate and subtropical zone, Henan is characterized by a humid and semi-humid monsoon climate featuring a cold winter, a dry spring, a hot summer and a clear autumn [10]. The annual average temperature is generally between 12 and 16 • C, and the annual precipitation varies from 530 to 1360 mm [11].
Winter wheat is one of main crops cultivated in Henan province, and is generally sown in October and harvested in May of the coming year. In this study, we are able to collect the data for yield per unit area of winter wheat (kg/ha) during the period 1988-2017 in 17 cities of Henan Province (Figure 2) from the Agricultural Statistical Yearbooks. The data for Jiyuan city were excluded in this study because of the limited period of available data (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). At the scale of whole province, the records of several indicators for agricultural productivity (APIs) and disasters statistics (DSIs) (listed in Table    , indicating that planting of winter wheat in these regions was characterized by higher average reduction rate and higher frequency of serious and extreme reduction years. According to the gradient of PCS1, Sanmenxia and Luoyang were classified as the highest-risk areas of CIYL, followed by Shangqiu, Zhoukou and Xinyang, then Zhenzhou, Luohe and Zhumadian. Similarly, the panel nonlinear regression analysis between wheat yields and climatic variables indicated that wheat planting in westernmost Henan was more sensitive to changes in the climate [19]. The westernmost areas are character-

Theory of Climate-Induced Yield
Following the theory of yield separation, the climate-induced yield (CIY) of crops can be calculated as: where AY i and TY i are the actual yield and trend yield of ith year, respectively, and TY are the fitting results by applying commonly-used simulation methods (e.g., 5-year moving average, linear regression, HP filter, etc.) on the AY series [11]. CIY usually varies with the fluctuation of climatic conditions, and TY is mainly determined by the development level of productive forces, such as contemporaneous technology, policy, and any other non-climatic factors [12]. Due to the similar agricultural development processes, the simulation method for TY within a certain province should generally be the same. Based on the CIY series, negative CIY values were considered as absolute and positive CIY values were assigned as zero, then the annual CIYL series were constructed. Seven Risk assessment indicators of yield losses (RAIYL) concerning the intensity, frequency, and variability of reduction were then extracted from the CIYL series; these are listed in Table 2.

Spatio-Temporal Analysis of CIYL and RAIYL
Principal component analysis (PCA) is a multivariate dimension reduction technique for turning a set of inter-correlated variables into new uncorrelated variables [13]. The new generated variables are often called principal components (PC), the linear combination of the original variables. The normalized eigenvectors (loadings) resulting from PCA represent the correlation between the original data and the corresponding PC. PCA can be further defined using six calculation modes based on the configuration of input data matrix; R-mode and S-mode PCA are most commonly applied to climatic series [14]. R-mode PCA (a data matrix with rows for the samples and columns for the related indicators) can simplify the interrelationship among the considered indicators, and is often used for comprehensive evaluation with multiple indicators [15]. S-mode PCA (a data matrix with rows for the time series and columns for the samples) can identify the homologous temporal patterns of a certain meteorological indicator in space [16]. In addition, the Mann-Kendall (M-K) test and Pettitt's test were applied in this study to recognize the changing features of annual CIYL series. The M-K test is a non-parametric trend test that is applicable for analyzing non-normally distributed data when the sample data are not necessarily compliant with a specific distribution [17]. Pettitt's test is widely used in identifying an abrupt change point in a time series and judging whether there is a significant difference in the cumulative distribution function before and after that point [18].

Calculations of CIYL for Winter Wheat
The crop TY series calculated by different yield decoupling methods can be expected to have different temporal characteristics. Hence, the selection of TY fitting methods has an important impact on the accuracy and objectivity of the subsequent calculation of CIY. The annual API series at the provincial level can intuitively reflect the inter-annual changes in productive forces, which is helpful for assessing the decoupling results. In this study, the wheat TY generated in each city by different fitting models were averaged, then the correlations between the various provincial TYs and the five APIs were calculated, as shown in Table 3. In view of the most significant correlations with APIs, the provincial TY generated using the quadratic curve model was more consistent with productive force changes in Henan Province. Therefore, the quadratic curve model was more suitable for extracting wheat TY series in Henan. On this basis, the annual CIYL series and RAIYL group in each city were then constructed. Table 3. Correlations between regional average TY generated by different fitting models and agricultural productivity indicators (** indicate the significant level at p < 0.01).

Spatial Pattern in Risk of CIYL
The statistical results of the seven RAIYL in each city are listed in Table 4. The average NYYR calculated for Henan province was 13.94, with the lowest coefficient of variation (CV) at 12.29%. In contrast, the values of CV for AYRR, NYYR_mo, NYYR_se, and NYYR_ex were all greater than 50%, especially for the CV of NYYR_se and NYYR_ex, which exceeded 100%. This may indicate that there were obvious spatial variations for RAIYL in Henan province. R-mode PCA was applied on the 17 × 7 matrix (17 refers to the number of cities, and 7 denotes the types of RAIYL). The first two PC explained 69.69% of the total variance, and the eigenvalues were higher than 1 (Table 5). Therefore, PC1 and PC2 can be effectively extracted in order to characterize the original group of indicators. According to the higher absolute values of the PC loadings in Table 5, PC1 mainly denotes the five indicators NYYR, AYRR, NYSYR_sl, NYYR_se, and NYYR_ex, while PC2 mainly denotes the CVC. The principal component scores (PCS) of each city shown in Figure 3 visually depict the spatial pattern of the risk of CIYL. Figure 3a shows that the positive values of PCS1 were Water 2021, 13, 3341 7 of 13 mainly found in the cities located in the westernmost, eastern, and southern areas of Henan province, indicating that planting of winter wheat in these regions was characterized by higher average reduction rate and higher frequency of serious and extreme reduction years. According to the gradient of PCS1, Sanmenxia and Luoyang were classified as the highest-risk areas of CIYL, followed by Shangqiu, Zhoukou and Xinyang, then Zhenzhou, Luohe and Zhumadian. Similarly, the panel nonlinear regression analysis between wheat yields and climatic variables indicated that wheat planting in westernmost Henan was more sensitive to changes in the climate [19]. The westernmost areas are characterized by rolling mountains, poor soil fertility, and a lower level of production technology, which may lead to the higher vulnerability of crop production to climate anomalies [20]. As shown in Figure 3b, the spatial pattern of PCS2 indicated that the CIYL in Zhumadian, and Xinyang had more obvious inter-annual variations.

Regional Changing Features of CIYL
S-mode PCA was applied on the 30 × 17 matrix (30 refers to the length of CIYL, and 17 denotes the number of cities). The first four PC explained 77.13% of the total variance, and the eigenvalues were larger than 1 (Figure 4). Therefore, the annual CIYL series during the period 1988-2017 in 17 cities can be simplified as four temporal modes represented by the mean of the cities of each cluster. The higher R 2 for the PCS-CIYL relationships in Figure 6 verifies the reliability of regional CIYLI for the subsequent studies.

Regional Changing Features of CIYL
S-mode PCA was applied on the 30 × 17 matrix (30 refers to the length of CIYL, and 17 denotes the number of cities). The first four PC explained 77.13% of the total variance, and the eigenvalues were larger than 1 (Figure 4). Therefore, the annual CIYL series during the period 1988-2017 in 17 cities can be simplified as four temporal modes represented by the first four PC. The higher loadings (>0.6) in Figure 5a make possible to divide Henan into four sub-regions with different CIYL variations. In Figure 5b, five cities in west Henan, five cities in north Henan, four cities in east Henan, and three cities in south Henan were classified as Regions I, II, III and IV, respectively. Similar geographical divisions were identified for the inter-annual variability of annual precipitation, mean temperature and minimum temperature across Henan Province [9,11]. These were mainly attributed to the multiple transitions within Henan: from south to north, there is a transition from humid to semi-humid monsoon climate, and from west to east, a transition from mountainous areas to plains areas [21,22]. Then, the regional CIYL were calculated from the arithmetic mean of the cities of each cluster. The higher R 2 for the PCS-CIYL relationships in Figure 6 verifies the reliability of regional CIYLI for the subsequent studies.
S-mode PCA was applied on the 30 × 17 matrix (30 refers to the length of CIYL, and 17 denotes the number of cities). The first four PC explained 77.13% of the total variance, and the eigenvalues were larger than 1 (Figure 4). Therefore, the annual CIYL series during the period 1988-2017 in 17 cities can be simplified as four temporal modes represented by the mean of the cities of each cluster. The higher R 2 for the PCS-CIYL relationships in Figure 6 verifies the reliability of regional CIYLI for the subsequent studies.   The correlations between the CIYL of each sub-region and provincial DSIs are shown in Table 6. These correlations can identify the leading meteorological disasters of winter wheat. In Region I, significant correlations were observed between CIYL and DD1/DD2, which may show that drought was the main meteorological disaster causing yield losses. As for Regions II and III, the significant correlations between CIYL and FD1/FD2 indicate that freezing had the most significant impact on yield fluctuations. In contrast, CIYL in Region IV were more significantly correlated to WD1 and WD2, which shows that yield losses were mainly caused by waterlogging. All of these results are in agreement with the geographical divisions for flood/drought and low-temperature disaster risk in Henan, which mainly result from the southeast-to-northwest decreasing gradient of precipitation and the south-to-north decreasing gradient of temperature [22].  The correlations between the CIYL of each sub-region and provincial DSIs are shown in Table 6. These correlations can identify the leading meteorological disasters of winter wheat. In Region I, significant correlations were observed between CIYL and DD1/DD2, which may show that drought was the main meteorological disaster causing yield losses. As for Regions II and III, the significant correlations between CIYL and Figure 6. Annual PCS series and regional average CIYL corresponding to Region I (a), Region II (b), Region III (c), Region IV (d) during study period (R 2 denotes the coefficient of determination for linear relationship between CIYLI and PCS). Table 6. Correlations between CIYL in each sub-region and provincial DSI (** and * indicate the significant level at p < 0.01 and 0.05, respectively).  Figure 7 shows the detection results for annual CIYL series in each sub-region. All the CIYL exhibited non-significant decreasing trends, indicating a weakened tendency toward risk of yield losses during the past three decades. According to the M-K statistics value, the decreasing trends of CIYL in Regions I and II were more notable. Furthermore, Figure 7 also shows the abrupt changes in the CIYL series which were detected using Pettitt's test. This shows that a significant declining trend occurred since the 2000s in each sub-region. The point of abrupt change in Regions I, II, III and IV was in 2002, 2006, 2005, and 2005, respectively, and all the detection results in Region I, II, and IV passed the significance test. These results are consistent with the findings of Shi and Tao [2], who indicated that the impacts of climate disasters in the main wheat production provinces of north China decreased over a range of decades.

Responses of CIYL Fluctuations to LOACI
In each sub-region, the correlation coefficients between CIYL and LOACI in different months were calculated (Table 7). In the months of the non-growing period of winter wheat, the CIYL-LOACI correlations ranged from −0.463 to 0.464. The TNA of July (TNA_Jul_cu), TNA of May (TNA_May_cu), Nino 1+2 of September (Nino 1+2_Sep_cu), and PNA of August (PNA_Aug_cu) in the current year had the strongest influence on CIYLI for Regions I, II, III, and IV, respectively. Concerning the months of the growing period, the CIYL-LOACI correlations ranged from −0.400 to 0.625. The EP/NP of March (EP/NP_Mar_co), SOI of February (SOI_Feb_co), WHWP of February (WHWP_Feb_co), and SOI of May (SOI_May_co) in the coming year had the strongest influence on CIYL for Regions I, II, III, and IV, respectively. It is noteworthy that the correlation between CIYL and WHWP_Feb_co in east Henan (Region III) was most significant, with a maximum value of 0.625. In addition, the R 2 for linear relationship between CIYL and two key LOACI in east Henan was higher when compared to those in other sub-regions. West, north, and south Henan are characterized by complex mountainous and hilly topography, which may lead to less predictable climate anomalies.
Water 2021, 13, x FOR PEER REVIEW 5 of 5 FD1/FD2 indicate that freezing had the most significant impact on yield fluctuations. In contrast, CIYL in Region IV were more significantly correlated to WD1 and WD2, which shows that yield losses were mainly caused by waterlogging. All of these results are in agreement with the geographical divisions for flood/drought and low-temperature disaster risk in of climate disasters in the main wheat production provinces of north China decreased over a range of decades.   Table 7. Summary of the most significant correlations between CIYL and LOACI during different months in each sub-region (gray fill color denotes the types of key LOACI before and during the growing period; _cu and _co denote the current year and coming year, respectively).

Region I Region II Region III Region IV
The current year Wheat yield can be considered as the sum of trend yield and climate-induced yield. In this study, TY was significantly correlated with year, and CIY was significantly correlated with LOACI. Hence, a multiple linear regression model including year and LOACI should provide an effective and convenient prediction method for actual yield. The simulation results were assessed by using Kling-Gupta efficiency (KGE), root-mean-square error (RMSE), and average relative error (ARE). KGE is an integrated evaluation index utilizing c, α, and β where c is the correlation coefficient between observed and simulated records, α is the standard deviation of simulated records over the standard deviation of observed records, and β is equal to the mean of simulated records over the mean of observed records [23]. As displayed in Table 8, the values of KGE for a single-factor model with predictors only years ranged from 0.801 to 0.915, while the KGE for multifactorial models with predictors of both years and various LOACI were higher (ranging from 0.882 to 0.958). These results indicate that the application of LOACI as an additional input greatly improved yield prediction. Among the four sub-regions, the highest prediction accuracy of the multifactorial model was found in Region II, with an ARE of only 3.484%, while the lowest was found in Region IV, with an ARE of 6.874%.

Discussion
Since the 2000s, significant reductions in climate-induced yield losses for winter wheat were found across Henan province. Similarly, panel regression analysis at the provincelevel in the Huang-Huai-Hai Plain showed that rising temperature resulted in increased yield in Henan province [24]. The significant increase of average and minimum temperature in the wheat growing season which occurred after 2000 may have reduced the frequency of low temperature disasters [19]. In addition, the number of rainy and flood days during spring in the Huang-Huai-Hai Plain showed a long-term decrease, especially after 1998, reducing the negative influence of excess precipitation on yields [6].
Previous studies of climate-induced yield of winter wheat in Henan indicate that the correlations between CIY and various agro-meteorological indicators ranged from −0.39 to 0.32, and the optimum starting months for CIY prediction were set in spring [25,26]. In contrast, the relationships between CIYL and LOACI in this study were more significant, with correlations ranging from −0.34 to 0.63. In addition, LOACI in both the growing period and non-growing period were significantly correlated to CIYL, which can contribute to the structuring of an early warning system of loss risk with flexible starting times. Especially, such outliers of antecedent signals as TNA_Jul_cu, TNA_May_cu, Nino 1 + 2_Sep_cu, and PNA_Aug_cu would seem to indicate a higher loss risk. These may be mainly attributed to the significant impact of LOACI on winter-spring climatic variations, with a lag time of 6-12 moths [27]. For example, a higher PNA in summer would indicate excess spring rainfalls in north China during the coming year, while a higher Nino 1 + 2 in Autumn would indicate a higher frequency of low-temperature extremes in northern China during the coming year [28,29]. Furthermore, linear regression with the independent variables of year and various LOACI proved to be a fast and simple forecasting method of actual wheat yields. However, crop yield variability is also associated with technology, genetics, climate, soil, field management, fertilizer application, [9] etc. In future works, it should be possible to achieve higher prediction accuracy in the yield forecasting method by including agricultural big data, circulation signals, and machine learning.

Conclusions
The results from applying R-mode PCA to the seven risk assessment indicators of yield losses (RAIYL) in 17 cities of Henan province indicate that the planting of winter wheat in the westernmost and eastern areas had a high risk of climate-induced yield losses (CIYL), especially in the cities of Sanmenxia and Luoyang. Using R-mode PCA on the CIYL series in 17 cities, four temporal modes from 1988-2017 were extracted; the corresponding sub-regions were west, north, east, and south Henan province. TNA, Nino 1 + 2, PNA, EP/NP, SOI, and WHWP were identified as the key LOACI affecting CIYLI. Compared to the other three sub-regions, the CIYL-LOACI relationship in east Henan was the most significant; a higher Nino 1 + 2 in September of the current year and a higher WHWP in February of the coming year can reliably indicate yield reductions in east Henan. Furthermore, linking the yield with LOACI and year could provide a convenient and efficient yield forecasting method.
Author Contributions: H.Z. was mainly committed to data analysis/interpretation and manuscript preparation; J.H. was mainly committed to study design and review; J.C. was mainly committed to data acquisition. All authors have read and agreed to the published version of the manuscript.