The Different Spatial and Temporal Variability of Terrestrial Water Storage in Major Grain-Producing Regions of China

: Irrigation is an important factor affecting the change of terrestrial water storage (TWS), especially in grain-producing areas. The Northeast China Plain (NECP), the Huang-Huai-Hai Plain (HHH) and the middle and lower reaches of the Yangtze River Basin Plain (YRB) are major grain-producing regions of China, with particular climate conditions, crops and irrigation schemes. However, there are few papers focusing on the different variation pattern of water storage between NECP, HHH and YRB. In this paper, the characteristics of terrestrial water storage anomaly (TWSA) and groundwater storage in the three regions mentioned above from 2003 to 2014 were analyzed, and the main reasons for water storage variations in the three regions were also discussed. The result shows that although effective irrigated areas increased in all three regions, TWSA only decreased in HHH and TWSA in the other two regions have shown an increasing trend. Spatially, the water storage deﬁcit was more serious in middle and south NECP and HHH. In the three regions, water storage variations were impacted by meteorological condition and anthropogenic stress (e.g., irrigation). However, irrigation water consumption has a greater impact on water storage deﬁcit in HHH than the other two regions, and water storage variation in YRB was mainly impacted by meteorological conditions. In this case, we suggest that the structure of agricultural planting in HHH should be adjusted to reduce the water consumption for irrigation.


Introduction
Freshwater is a significant resource to sustain agriculture and ensure food security. Water consumption in irrigation accounts for about 70% of global freshwater [1,2]. Terrestrial water storage (TWS), as an important component of water cycle, is strongly impacted by irrigation. Because of the diversity of climate conditions and human activity intensity, TWS changes are different from various regions [3]. The Northeast China Plain, the Huang-Huai-Hai Plain HHH and the middle-lower reaches of the Yangtze River Basin (YRB) are three major grain-producing regions of China with diverse climate conditions and agriculture planting characteristics. Thus, what is the difference of TWS variability in the three grain-producing areas?
TWS is defined as all forms of water stored above and underneath the surface of the Earth, and is a combination of snow, ice, soil moisture (SM), groundwater, surface water and biomass [4,5]. At present, TWS change is monitored in three ways, including in situ observations, global hydrological models and satellite remote sensing. Because of

Study Area
The Northeast China Plain (NECP) is located in Northeast China (Figure 1a), with an area of 728,000 km 2 (and the crop land accounts for 38.7%). It belongs to the temperate monsoon climate, which experiences hot and humid conditions in summer and dry and cold conditions in winter. The annual precipitation ranges from 350 to 1200 mm and the mean annual temperature ranges from −4.8 to 11.3 • C [28]. The cropping system brings in one harvest per year and its main agricultural products are rice, maize and soybean. The Huang-Huai-Hai Plain (HHH) is located in North China (Figure 1b), with an area of 508,000 km 2 ; the crop land accounts for 59%. The mean annual precipitation ranges from 600 to 800 mm and the annual accumulated temperature is about 4500 • C [29]. The main crops are wheat and maize, which can be harvested three times every two years. The middle and lower reaches of the Yangtze River Basin Plain (YRB) (Figure 1c) is located in a subtropical zone with an area of 916,000 km 2 ; the crop land accounts for 37.8%. The annual precipitation ranges from 1100 to 1500 mm. The main crop is rice and the cropping system is yield two crops a year. Crops in NECP and HHH are mostly fed by groundwater, which has resulted in a groundwater deficit. Meanwhile, TWS changes in YRB is influenced both by climate change and human activities (i.e., the Three Gorges Dam, TGD) [30,31].

Datasets
The GRACE dataset was provided by the Center for Space Research at the University of Texas at Austin (CSR). Spherical harmonics for CSR RL05 were truncated at the maximum degree and order of 60, detriped and filtered using a 300 km Gaussian filter method to suppress GRACE measurement noise of high-degree and order spherical harmonics [32]. Gridded TWS anomalies (TWSA) data, which were derived from GRACE data CSR RL05 and resampled at a spatial resolution of 1 × 1 • , were used to evaluate water storage deficit (WSD), covering the period from January 2003 to December 2014. The GRACE mission was the first to allow the monitoring of large-scale TWS anomalies [33]. Using this new information, in combination with external datasets, many studies analyzed the temporal evolution of GWS anomalies at regional and basin scales [34,35]. What needs to be explained is that the spatial resolution of GRACE data are coarse, and detailed spatial information of water storage change cannot be presented. The self-calibrating Palmer drought severity index (scPDSI) dataset and monthly precipitation dataset were used to characterize climate conditions of the three areas. The scPDSI dataset was developed by Wells et al. [36] and provided by the Climatic Research Unit at the University of East Anglia (CRU) at a spatial resolution of 0.5 × 0.5 • . The weight coefficient and duration factor of scPDSI for each location is estimated based on the observations of that location only and it is more appropriate for geographical comparison [37]. The scPDSI data are appropriate for China according to [38]. The monthly precipitation dataset was provided by the National Meteorological Information Center (NMIC) of the Chinese Meteorological Administration (CMA). It is developed by the Thin Plate Spline (TPS) method based on 2472 meteorological stations in China at a spatial resolution of 0.5 × 0.5 • .

Datasets
The GRACE dataset was provided by the Center for Space Research at the University of Texas at Austin (CSR). Spherical harmonics for CSR RL05 were truncated at the maximum degree and order of 60, detriped and filtered using a 300 km Gaussian filter method to suppress GRACE measurement noise of high-degree and order spherical harmonics [32]. Gridded TWS anomalies (TWSA) data, which were derived from GRACE data CSR RL05 and resampled at a spatial resolution of 1 × 1°, were used to evaluate water storage deficit (WSD), covering the period from January 2003 to December 2014. The GRACE mission was the first to allow the monitoring of large-scale TWS anomalies [33]. Using this new information, in combination with external datasets, many studies analyzed the tem- The yearly irrigation water consumption (IWC) data and monthly GWS data, which were calculated based on WaterGAP 2.2 at a spatial resolution of 0.5 × 0.5 • , were used to evaluate the impact of irrigation and groundwater storage variation, respectively. Here, water consumption is defined as the part of withdrawn water that does not return to the system but is rather evaporated or incorporated in products [8]. Temporal GWS changes are calculated as follows: where dS GW dt is the groundwater storage changed over time t, R g is the groundwater recharge, Q g is the base flow and NA g is the net extraction of the groundwater.
Information of datasets used in this study are shown as Table 1. Note that here, the "-" represents that the dataset was from the authors Müller et al. [8,39] and Döll et al. [40].

Methods
TWSA, derived from GRACE data, is used to calculate the water storage deficit and detect the difference in TWS variability. The interannual trends were derived from Seasonal and Trend decomposition using Loess (STL).

Water Storage Deficit Index
The monthly water storage deficit (WSD), which is defined as the difference between TWSA in any month and the average values, is used to measure hydrological drought characterization [19,25], as per Equation (2): where WSD i,j represents the WSD for the jth month in year i, TWSA i,j is the GRACEinferred TWSA time series for the jth month in year i and TWSA j is the long-term mean (from January 2003 to December 2014) of TWSA for the same month. Negative WSD represents that TWS is less than the normal value. To better characterize droughts based on WSD, Sun et al. [25] normalized WSD into the WSD index (WSDI) by the zero mean normalization method as Equation (3): where µ and σ are the mean and standard deviation of the WSD time series, respectively. The magnitude of WSDI indicates the drought intensity [25].

Groundwater Storage Anomaly
GWSA is defined as the difference between GWS in any month and the average value. GWSA is used to indicate groundwater storage variation characterization, as per Equation (4): where GWS i,j represents the GWS for the jth month in year i and GWS is the long-term mean (from January 2003 to December 2014) of GWS for the same month.

Time Series Decomposition
The time series was composed of seasonal signal (S seasonal ), trend signal (S trend ) and residuals based on the non-parametric STL approach (as per Equation (5)) [41]. STL decomposition have been applied to GRACE data to identify the signatures of extreme climate events such as droughts and floods on TWS [25,42].
where S total represents the original signal. The STL method is a robust and computationally efficient approach commonly used for detecting non-linear patterns in trend estimates [25]. In this work, STL approach is used for reducing the influence of seasonal variation on long-term trends. Interannual trends of TWSA and WSDI were extracted by performing the 'data_decomp' function provided by R and long-term trends were estimated by the linear regression method.

Validation of WSDI
WSDI is a reliable indicator of drought [25]. WSDI shows changes in water storage more intuitively. In this paper, scPDSI was used as a reference to validate WSDI derived from GRACE production. Figure 2 shows the comparison of scPDSI and WSDI from 2003 to 2014 in the three major grain-producing regions. The variability of WSDI agreed well with that of scPDSI with correlation coefficients of 0.71, 0.69 and 0.41, respectively. However, there were also obvious differences between variations derived from scPDSI and WSDI. For instance, in NECP and HHH, WSDI was lower than scPDSI during the wet period, while it was higher than scPDSI during the drought period. In YRB, the characteristics of WSDI was obviously different from that of scPDSI in 2003, even though they agreed well with each other during 2004-2014, with a correlation coefficient of 0.73.
Possible reasons for those differences are as follows: (1) WSDI and scPDSI are calculated from different data sources and through different methods. For instance, WSDI is derived from TWSA, and scPDSI is estimated based on meteorological datasets; (2) resampling is required during the analysis, since the spatial resolutions of the two measurements are different (e.g., the spatial resolution of WSDI is 1 × 1 • , while that of scPDSI is 0.5 × 0.5 • ); (3) WSDI is impacted not only by meteorological conditions but also by anthropogenic factors, while scPDSI is impacted only by meteorological conditions.

The Temporal Variation of Water Storage in Major Grain-Producing Regions
The interannual trend of TWSA was derived by STL method, while the seasonal signal was removed to reduce the impact to long-term trend. Figure  Moreover, TWSA has changed the most in HHH during the past 12 years among the three grain-producing regions. Figure 4 shows the annual variation of TWSAs in the three regions, and it can be seen that the annual differences in the three regions were relatively obvious. TWSA shows two crests in March and August, and two troughs in June and October in NECP, while TWSA minimized in June in HHH and peaked in July in YRB. In addition, compared with the drastic changes of TWSA in both YBR and HHH, there was no significant intra-annual change in NECP. It is worth noting that TWSA in HHH have shown a loss throughout the year.

Water Storage Deficit in Major Grain-Producing Regions
WSD is an important characteristic of drought. Figure  Possible reasons for those differences are as follows: (1) WSDI and scPDSI are ca lated from different data sources and through different methods. For instance, WSD derived from TWSA, and scPDSI is estimated based on meteorological datasets; resampling is required during the analysis, since the spatial resolutions of the two m urements are different (e.g., the spatial resolution of WSDI is 1 × 1°, while that of scP  Figure 4 shows the annual variation of TWSAs in the three regions, and it c that the annual differences in the three regions were relatively obvious. TWSA crests in March and August, and two troughs in June and October in NECP, w minimized in June in HHH and peaked in July in YRB. In addition, compare drastic changes of TWSA in both YBR and HHH, there was no significant in change in NECP. It is worth noting that TWSA in HHH have shown a loss thro year.    Figure 4 shows the annual variation of TWSAs in the three regions, and it can be seen that the annual differences in the three regions were relatively obvious. TWSA shows two crests in March and August, and two troughs in June and October in NECP, while TWSA minimized in June in HHH and peaked in July in YRB. In addition, compared with the drastic changes of TWSA in both YBR and HHH, there was no significant intra-annual change in NECP. It is worth noting that TWSA in HHH have shown a loss throughout the year.     Cumulative WSD indicates the continuity of the water storage deficit. Accor Sun et al. [25], a declining trend of cumulative WSD represents a lasting deficit, w increasing trend represents a surplus. In NECP (as shown in Figure 5a Figure 5c). Furthermore, these points of WSD in NECP, HHH and YRB are consistent with that of the WSD in th regions accordingly. Figure 6 presents time variations in WSDI in the period of 2003-2014 and the s trends were removed from WSDI time series by STL. As can be seen, the trend o was consistent with that of TWSA (shown in Figure 3) in three grain-producing accordingly. WSDI in NECP increased slightly from 2003 to 2014, with a rate of shown in Figure 6a Cumulative WSD indicates the continuity of the water storage deficit. According to Sun et al. [25], a declining trend of cumulative WSD represents a lasting deficit, while an increasing trend represents a surplus. In NECP (as shown in Figure 5a), increasing trends were detected during 2003-2007 and 2013-2014, while a declining trend was found during 2007-2012. In addition, the cumulative WSD in NECP has been negative since 2008, which indicated the continued loss of water storage since then. As shown in Figure 5b, the cumulative WSD in HHH was positive during 2004-2014, but it has declined since August 2008. The cumulative WSD in YRB was generally negative during the whole period, which was consistent with the TWSA. The cumulative WSD in YRB decreased from 2003 to July 2008, and then increased until 2014 (as shown in Figure 5c). Furthermore, these turning points of WSD in NECP, HHH and YRB are consistent with that of the WSD in the three regions accordingly. Figure 6 presents time variations in WSDI in the period of 2003-2014 and the seasonal trends were removed from WSDI time series by STL. As can be seen, the trend of WSDI was consistent with that of TWSA (shown in Figure 3) in three grain-producing regions accordingly. WSDI in NECP increased slightly from 2003 to 2014, with a rate of 0.03 (as shown in Figure 6a Figure 6b,c, WSDI declined with a rate of 0.17 and increased with a rate of 0.08, respectively. The turning points of WSDI in HHH (from drought to wet) and YRB (from wet to drought) appeared during 2008-2009, which is consistent with that of WSD, shown in Figure 5.  Fig  WSDI declined with a rate of 0.17 and increased with a rate of 0.08, respectively. ing points of WSDI in HHH (from drought to wet) and YRB (from wet to dro peared during 2008-2009, which is consistent with that of WSD, shown in Figur

The Spatial Pattern of WSDI Variation in Major Grain-Producing Regions
Annual average of WSDI from 2003 to 2014 were also calculated in order t the spatial distribution of WSDI variation in major grain-producing regions. A seen from Figure

The Spatial Pattern of WSDI Variation in Major Grain-Producing Regions
Annual average of WSDI from 2003 to 2014 were also calculated in order to analyze the spatial distribution of WSDI variation in major grain-producing regions. As can be seen from Figure  As shown in Figure 7, WSDIs were negative in north NECP; however, it can be seen from Figure 6a that WSDI was positive in 2005. It suggested that the water storage deficit was not significant in north NECP in 2005. Furthermore, water storage deficit was serious in middle NECP in 2009 and 2011, as shown in Figure 6, although the water storage was surplus (as shown in Figure 5a). In the junctional zone of HHH and YRB, water storage deficit variations were different from other regions of YRB, but consistent with water storage deficit variations in HHH. The possible reasons for this difference are as follows: (1) the junctional region is located at the same latitude as south HHH and the climate conditions are similar; (2) the geographic elements (i.e., terrain, river system) in the junctional region are similar to south HHH; (3) water storage variations may be impacted by spatial autocorrelation. As shown in Figure 7, WSDIs were negative in north NECP; however, it can be seen from Figure 6a that WSDI was positive in 2005. It suggested that the water storage deficit was not significant in north NECP in 2005. Furthermore, water storage deficit was serious in middle NECP in 2009 and 2011, as shown in Figure 6, although the water storage was surplus (as shown in Figure 5a). In the junctional zone of HHH and YRB, water storage deficit variations were different from other regions of YRB, but consistent with water storage deficit variations in HHH. The possible reasons for this difference are as follows: (1) the junctional region is located at the same latitude as south HHH and the climate conditions are similar; (2) the geographic elements (i.e., terrain, river system) in the junctional region are similar to south HHH; (3) water storage variations may be impacted by spatial autocorrelation.

GWS variations in Major Grain-Producing Regions
As mentioned above, the groundwater storage is an important component of TWS.

GWS Variations in Major Grain-Producing Regions
As mentioned above, the groundwater storage is an important component of TWS. Figure 8 shows that GWS was stable in NECP and YRB, while it decreased significantly in HHH, during 2003-2014. According to the average of GWSA (shows as Table 2), GWS was in serious deficit in HHH, with an average of −194.70 cm and a minimum of −288.64 cm. In contrast, GWS in YRB was relatively abundant, with an average of 72.94 cm and a minimum of 41.77 cm. The GWS in NECP was in the middle of the three grain-producing regions, with an average of 5.03 cm.
In addition, the maximum and minimum GWS in NECP and YRB correspond well to floods and drought events, such as the flood disaster in the Songhua River Basin in August 2013 and the severe drought events in the middle and lower reaches of the Yangtze River in 2003-2004. It can be seen from the spatial distribution of GWS change trend (as shown in Figure 9) that GWS in the northern part of NECP was increasing, while in the southern and central parts it showed a decreasing trend. The decreasing area accounted for about 35% of the total area of NECP. GWS had an increasing trend only in the northeast of HHH, and the decreasing area accounted for about 88% of the total area of HHH. GWS in the eastern part of YRB increased significantly, and the reduction areas were mainly distributed in the junction area with HHH; the reduction area accounted for about 22% of YRB.  In addition, the maximum and minimum GWS in NECP and YRB correspond to floods and drought events, such as the flood disaster in the Songhua River Bas August 2013 and the severe drought events in the middle and lower reaches of the Y tze River in 2003-2004. It can be seen from the spatial distribution of GWS change t (as shown in Figure 9) that GWS in the northern part of NECP was increasing, wh the southern and central parts it showed a decreasing trend. The decreasing are counted for about 35% of the total area of NECP. GWS had an increasing trend only i northeast of HHH, and the decreasing area accounted for about 88% of the total ar HHH. GWS in the eastern part of YRB increased significantly, and the reduction were mainly distributed in the junction area with HHH; the reduction area accounte about 22% of YRB.    In addition, the maximum and minimum GWS in NECP and YRB correspond well to floods and drought events, such as the flood disaster in the Songhua River Basin in August 2013 and the severe drought events in the middle and lower reaches of the Yangtze River in 2003-2004. It can be seen from the spatial distribution of GWS change trend (as shown in Figure 9) that GWS in the northern part of NECP was increasing, while in the southern and central parts it showed a decreasing trend. The decreasing area accounted for about 35% of the total area of NECP. GWS had an increasing trend only in the northeast of HHH, and the decreasing area accounted for about 88% of the total area of HHH. GWS in the eastern part of YRB increased significantly, and the reduction areas were mainly distributed in the junction area with HHH; the reduction area accounted for about 22% of YRB.

Discussion
According to the water balance principle, water storage change is simultaneously impacted by meteorological condition and anthropogenic stress. Figure 10 shows the relationships between precipitation, scPDSI and TWSA in the three major grain-producing regions. As shown in Figure 10, the meteorological condition tended toward wet in NECP

Discussion
According to the water balance principle, water storage change is simultaneously impacted by meteorological condition and anthropogenic stress. Figure 10 shows the relationships between precipitation, scPDSI and TWSA in the three major grain-producing regions. As shown in Figure 10, the meteorological condition tended toward wet in NECP and YRB (Figure 10a,c), while it tended to drought in HHH (Figure 10b). Moreover, it suggested that variations of precipitation were consistent with that of scPDSI and WSDI (according to Figure 10). The annual average precipitation fluctuated greatly after 2009 in NECP and YRB. In addition, annual variations of TWSA and precipitation were inconsistent in the three regions. On the one hand, TWSAs in YRB showed approximately one-month shorter lag to precipitation and TWSAs in HHH presented a time lag of two months to precipitation. while TWSAs showed two crests (in March and August) and two troughs (in June and October) yearly in NECP with only one peak in annual precipitation. This suggests that precipitation is not the most important reason for water storage change in NECP and HHH.
On the other hand, the seasonal variations of TWSA were different in the three regions (according to Figure 4). This may be caused by different agricultural irrigation and cropping systems; for example, in NECP, the main crops are ripe once in a year and the growing period is from April to October. Therefore, TWSA decreased from April to October, except for during the rainy season (from June to August) (Figure 10d). Similarly, TWSA decreased during the growing period of the main crops (winter wheat) in HHH but increased during the rainy season (Figure 10e). However, there was not significant decline in TWSA during the growing period of crops in YRB (Figure 10f), because of its abundant water resources. In general, variations of TWSA and precipitation were consistent with each other during 2003-2014, but the seasonal variations of them were different.
Furthermore, effective irrigation areas were increased in the three regions ( Figure 11). However, IWCs did not vary with the effective irrigation area. In fact, IWC varied with meteorological conditions. Figure 12 shows the relationship between precipitation and IWC in the three major grain-producing regions. As can be seen, IWC was reduced as compared to a normal year when precipitation was plentiful (annual precipitation anomaly was positive); conversely, IWC was increased to assure grain production when precipitation was poor (the value is negative). For example, HHH is the important producing region of winter wheat and irrigation is the necessary way to ensure grain output. IWC increased in HHH because of the lack of precipitation, leading to the water storage deficit getting more serious during 2003-2014, according to Figures 5, 6 and 8. This indicated that under the superposition of poor precipitation and increased IWC, water storage reduction and water storage deficit were further aggravated, especially in HHH [43]. Thus, we suggest that the structure of agricultural planting in HHH should be adjusted to reduce the water consumption for irrigation in winter.
Water 2021, 13, x FOR PEER REVIEW and IWC in the three major grain-producing regions. As can be seen, IWC was r as compared to a normal year when precipitation was plentiful (annual preci anomaly was positive); conversely, IWC was increased to assure grain productio precipitation was poor (the value is negative). For example, HHH is the importa ducing region of winter wheat and irrigation is the necessary way to ensure grain IWC increased in HHH because of the lack of precipitation, leading to the water deficit getting more serious during 2003-2014, according to Figure 5, 6 and 8. Th cated that under the superposition of poor precipitation and increased IWC, water reduction and water storage deficit were further aggravated, especially in HH Thus, we suggest that the structure of agricultural planting in HHH should be a to reduce the water consumption for irrigation in winter.

Conclusions
Due to differences in the meteorological condition and cropping systems, th storage change characteristics vary from region to region. Based on GRACE CSR d analyzed variation features of TWSA in three major grain-producing regions o over the period from 2003 to 2014. By calculating WSD and WSDI, we evaluate  and IWC in the three major grain-producing regions. As can be seen, IWC was reduced as compared to a normal year when precipitation was plentiful (annual precipitation anomaly was positive); conversely, IWC was increased to assure grain production when precipitation was poor (the value is negative). For example, HHH is the important producing region of winter wheat and irrigation is the necessary way to ensure grain output. IWC increased in HHH because of the lack of precipitation, leading to the water storage deficit getting more serious during 2003-2014, according to Figure 5, 6 and 8. This indicated that under the superposition of poor precipitation and increased IWC, water storage reduction and water storage deficit were further aggravated, especially in HHH [43]. Thus, we suggest that the structure of agricultural planting in HHH should be adjusted to reduce the water consumption for irrigation in winter.

Conclusions
Due to differences in the meteorological condition and cropping systems, the water storage change characteristics vary from region to region. Based on GRACE CSR data, we analyzed variation features of TWSA in three major grain-producing regions of China over the period from 2003 to 2014. By calculating WSD and WSDI, we evaluated water storage deficits both temporally and spatially. Then, we validated WSDI by using scPDSI production. Finally, we discussed the relationship between TWSA variations and indicators and the possible main reason that caused the different variation patterns. We found

Conclusions
Due to differences in the meteorological condition and cropping systems, the water storage change characteristics vary from region to region. Based on GRACE CSR data, we analyzed variation features of TWSA in three major grain-producing regions of China over the period from 2003 to 2014. By calculating WSD and WSDI, we evaluated water storage deficits both temporally and spatially. Then, we validated WSDI by using scPDSI production. Finally, we discussed the relationship between TWSA variations and indicators and the possible main reason that caused the different variation patterns. We found that TWSA increased in NECP and YRB with different rates but decreased in HHH during 2003-2014. Moreover, TWSA changed continuously in HHH and YRB while volatility was seen in NECP. In NECP, WSD had two obvious turning points in April 2007 and August 2012; in HHH, WSD has been continuous since 2008, while in YRB, the water storage has showed a surplus since 2008. Spatially, WSDI represented significant differences over NECP. WSD was more serious in middle and south NECP than in the north. The spatial heterogeneity of WSDI was in apparent in HHH and YRB. Furthermore, in the junctional region of YRB and HHH, WSD variations were consistent with that in HHH, because of the similar geography and climate conditions and the impact of adjacent areas. Water storage change is impacted both by meteorological condition and anthropogenic stress. In the three major grain-producing regions, variations in TWSA and precipitation show similar trends during 2003-2014, though their seasonal change characteristics were different. The water storage deficit in HHH was the most serious, since the IWC has been increasing. Last but not the least, increasing IWC to assure grain production under poor precipitation may cause serious WSD. In this case, we suggest that the structure of agricultural planting in HHH should be adjusted to reduce the water consumption for irrigation in winter.  Data Availability Statement: This study did not report any data.