Estimating the potential yield and ETc of winter wheat across Huang-Huai-Hai Plain in the future with the modified DSSAT model

The DSSAT model, integrated the calibrated Hargreaves ET model and dynamic crop coefficient, was run with the generated weather data by SDSM4.2 and CanESM2 to predict the potential yield and crop water requirement (ETC) of winter wheat in the Huang-Huai-Hai Plain in China under RCP4.5 and RCP8.5 scenarios. The results showed that the spatial distribution of potential yield in the future under RCP4.5 and RCP8.5 were similar, characterized by an increasing trend from the northwest inland to the southeast coast. The spatial distribution of ETC decreased gradually from the Shandong Peninsula to the surrounding area, and the minimum ETC was observed in the southern part of Huang-Huai-Hai Plain. The potential yield, ETC, and effective precipitation during winter wheat growing seasons might increase in the future under RCP4.5, while irrigation water requirements (IWR) would decrease. Under RCP8.5, the effective precipitation during the wheat growing seasons decreased first and then increased. However, the potential yield, ETC, and IWR of winter wheat increased first and then decreased. This study can provide some scientific evidence to mitigate the negative effects of climate change on agricultural production and water use in the Huang-Huai-Hai Plain.

The Huang-Huai-Hai Plain is an important base for high-quality wheat production in China, as well as food security of China. However, precipitation in this region is less than 150 mm, which can not meet the water requirement of about 450 mm for winter wheat 1 . Therefore, scientific irrigation management is the key to ensure stable and high yield of grain crops as well as high water use efficiency. The accurate estimation of crop water requirement is the key to irrigation scheduling. The IPCC fifth global climate change assessment report showed that global average temperature of 2016-2035 may rise by 0.3-0.7 °C and climate warming is expected to continue throughout the 21 st century compared with 1986-2005 2 . Climate change will affect crop growth season 3 , planting boundaries 4 , and cropping systems 5,6 . Relative studies showed that mean seasonal evapotranspiration (ET) of winter wheat in the North China Plain under well-watered conditions gradually increased from the 1980s to 2000s 7 and the reference ET (ET 0 ) will also increase in the future under representative concentration paths 8,9 . This study aims to explore the change of winter wheat water requirement (ET C ) under main climate scenarios in the future across Huang-Huai-Hai Plain. The results will provide a theoretical basis for the adjustment of agricultural production pattern, optimal allocation of water resources, and scientific response to climate change impact on agricultural production.
The estimation of ET C in the future has been taken more attention and the main works focus on formula and crop model methods. However, both methods need global climate models (GCMs) and downscaling methods to project local weather data. Formula method uses formula to estimate ET 0 and then crop coefficient is used to estimate ET C 10-12 . Due to climate change, crop growth period and physiological traits will change, and correspondingly, the value of crop coefficient will also change 7 . While mostly, the change of crop coefficient was not considered in previous studies. Though in some studies, the thermal time method was used to estimate phenology 1

Results
Parameter estimation and verification of DSSAT-CERES-Wheat. The results of genetic parameters estimation for the cultivar of 'Bainong 207' were showed in Table 1. A comparison of crop coefficient curves was showed in Fig. 1. The original K C varied daily between 1.0 and 1.1, which were greater than 1.0 during the whole growth period. The modified K C varied daily between 0.4 and 1.136, which were fitting well with observed ones (Fig. 1). The growth of winter wheat under full irrigation in 2015 and 2016 were simulated by three DSSAT modification plans (Table 2). When the genetic parameters, soil data, field management, and weather data were the same, the simulated flowering and maturity dates by the three modification plans were the same as the observed ones. However, the simulated maximum leaf area index (MLAI) by Plan 1 was the smallest, while the absolute relative error of MLAI by Plan 2 and 3 was smaller. The absolute relative error of biomass (at maturity) and yield changed very little among the three plans. The simulation results of LAI and biomass by three plans were showed in Fig. 2. Three curves of LAI looked as the same one, similar results could be found for biomass. The reason may be that only the module of evapotranspiration in DSSAT was modified, and the other modules had no change. The simulation results for both LAI and biomass were better in the middle stage of winter wheat growth and poor in the late stage of winter wheat growth. R 2 and RMSE between the simulated LAI and observed ones by three plans were always the same, R 2 was 0.878, and RMSE was 1.01. R 2 between the simulated biomass and observed ones by three plans always was 0.982, while RMSE between the simulated biomass and observed ones was 1705 kg ha −1 by Plan 1, 1712 kg ha −1 by Plan 2, 1713 kg ha −1 by Plan 3, respectively.
Results of ET 0 at the 7 stages in the 2015-2016 season estimated by PT, P-M, HS, H method were showed in Fig. 3. Compared with the P-M method, ET 0 estimated by PT method were lower at all stages (except for Stage 3), ET 0 estimated by H method were higher at all stages (except for Stage 2), ET 0 estimated by HS method were slightly lower in Stage 1-3 (overwintering and greening) and slightly higher in Stage 4-6 (jointing, heading and grain-filling). Taking the ET 0 estimated by P-M method as the standard, ET 0 estimated by HS method was the most close to P-M method in Stage 1, 4, 5, 6 (overwintering, Jointing, Heading, Grain-filling), in Stage 2 (overwintering) for H method, and in stage 3 (greening) for PT method. In the whole growth periods (Stage 7), ET 0 estimated by HS method was the most close to P-M method. Therefore, HS method had a better performance for calculating ET 0 compared to the P-M method in Huang-Huai-Hai Plain.  The simulated ET C and observed ET C at different stages in the 2015-2016 season were showed in Fig. 4. Compared with the observed ET C in different stages, the simulated ET C by Plan 1 had the great difference, indicating both the default PT method and the lack of crop coefficient may cause much errors. At stages from 1 to 3, the simulated ET C by Plan 1 were 12, 5.5 and 8.9 mm higher than the observed ones, respectively. However, at stages 4 to 6, the estimated ET C by Plan 1 were 9.7, 17.9, 11.2 mm lower than the observed ones. The reason may be that the Plan 1 had no crop coefficient. ET C in the early stages of winter wheat (stages 1-3) estimated by Plan 2 and 3 was similar and close to the observed data. However, at the late three stages (stages 4-6), ET C estimated by the Plan 3 was higher than those in plan 2, and much close to the observed data with the absolute difference of −4.1 to 2.0 mm. The reason might be that the calibrated Hargreaves method mainly improve the accuracy of ET C simulation at the later stages of winter wheat by integrating the crop coefficient. Considering the good performance for ET C estimation in the whole growth period, this study selected Plan 3 to estimate ET C in the future. Spatial distribution of potential yield of winter wheat under future scenarios. The modified DSSAT model with Plan 3 was used to simulate the growth of winter wheat without water stress under main climate scenarios in the future in the Huang-Huai-Hai Plain. The spatial patterns of potential yield and ET C under future scenarios are presented on Figs 5 and 6.
The spatial distributions of simulated potential yields in 2020s, 2050s, and 2080s under RCP4.5 in the Huang-Huai-Hai Plain were similar to those under RCP8.5. Yields in both scenarios increased gradually from the northwest inland to the southeast coast (Fig. 5). The area with a potential yield less than 10000 kg ha −1 a −1 under RCP4.5 was in Beijing-Tanggu region, the border between Hebei and Shandong Provinces (Fig. 5a-c). However, under RCP8.5, the area with yield less than 10000 kg ha −1 a −1 distributed in the east-central Hebei, north-central Henan, and northwestern Shandong in 2020s (Fig. 5d), then gradually decreased to the junctions of Hebei, Shandong and Henan Provinces in 2050s (Fig. 5e) and extended to Zhengzhou and Kaifeng in Henan Province in 2080s (Fig. 5f). The area with a potential yield more than 11500 kg ha −1 a −1 under RCP4.5 was near the coast of Jiangsu Province in 2020s and 2050s (Fig. 5a,b), then it expanded to the coast of Jiangsu and  Spatial distribution of ETC of winter wheat under future scenarios. The spatial distribution of ET C in 2020s, 2050s, and 2080s under RCP4.5were also similar to those under RCP8.5 ( Fig. 6), both decreased gradually from the Shandong Peninsula to the surrounding area, and the minimum ETc, generally less than 400 mm a −1 , was observed in the southeast of Huang-Huai-Hai Plain in Jiangsu Province and the border between Henan and Anhui Provinces. Under RCP4.5, the area with the simulated ET C of 440-460 mm a −1 accounted for about 35% of the total study area, mostly distributed in the central of Huang-Huai-Hai Plain ( Fig. 6a-c). However, under RCP8.5, the area with the simulated ET C of 460-480 mm a −1 was the largest (Fig. 6e-g). The area with the ETc more than 500 mm a −1 under RCP4.5 only distributed near Haiyang in Shandong Province (Fig. 6a-c), while it expanded to Weifang and Yiyuan countries, in the central part of Shandong Province under RCP8.5 ( Fig. 6e-g).

Changes of potential yield, ETC and IWR under future scenarios. The changes of main factors of win-
ter wheat in Huang-Huai-Hai Plain in different periods were showed in Table 3. The estimated potential yield of winter wheat in the three time periods showed an increasing trend under RCP4.5, however, it increased first and then reached to the highest in 2050s, after that it decreased in 2080s under RCP8.5. The change patterns of ET C were consistent with that of potential yield in the future. That's probably because of the moderate radiation forcing and greenhouse gas concentration (2050s and 2080s under RCP4.5, 2050s under RCP8.5) were beneficial to the potential yield and ET C of winter wheat, while the excessive radiation forcing and greenhouse gas concentration (2080s under RCP8.5) had an inhibitory effect. The effective precipitation during the growing period of winter wheat showed an increase trend from 195 to 208 mm a −1 under RCP4.5 and while it decreased first and then increased under RCP8.5. Irrigation water requirements (IWR) was different from potential yield, ETc and effective precipitation. It decreased slightly from 2020s to 2080s under RCP4.5, while it increased first and then decreased under RCP8.5.

Discussions
Although the downscaling method to solve the scale mismatch problem between GCMs and crop model is well done 3,5,[13][14][15] , the accuracy of simulated meteorological factors needs to be improved. If the meteorological data with low precision were input into the crop models, the error would accumulate for ET C. To study the change of ET C of winter wheat in the future, the main content is to study the impacts of meteorological factors on ET C. Thus, accurate meteorological data are needed. So far, the method based on temperature is the most appropriate to predict ET C in the future.
Climate change can affect the growth mechanism of winter wheat and ET C . On one hand, the increasing temperature can lead to a higher water vapor deficit and shorter growth period. The increasing concentration of CO 2 may reduce leaf stomatal conductance of winter wheat, which consequently reduces the ET C 32 . On the other hand, CO 2 fertilization and enhanced photosynthesis can increase ET C of winter wheat 2,13 . Consequently, the change of ET C is not remarkable 33 , which is confirmed in this study ( Table 3). The greenhouse gas concentration under RCP4.5 is lower than the target level and tend to be stable, while it will increasing all the time under RCP8.5 34 . The ET C of winter wheat under RCP4.5 was higher than that under RCP8.5, indicated that the higher radiation forcing and greenhouse gas concentration under RCP8.5 must inhibit the ET C of winter wheat. In addition, ET 0 in Huang-Huai-Hai Plain in 21st century under RCP4.5 and RCP8.5 showed an increasing trend 7 . The trend of ET 0 is different from the trend of ET C , which may infer the complicated influence of climate change on growth period and crop coefficient of winter wheat.
In this study, the higher yield was distributed in Jiangsu and Shandong Peninsula, and the higher ET C was distributed in the Shandong Peninsula in the future. However, this trend was different from that of Mo et al. 32 , who used the VIP model to simulate yield and ET C of winter wheat in Huang-Huai-Hai region. Mo et al. found that   relatively higher yield distributed in the Henan Province from 2000 to 2008, but ET C of winter wheat in Shandong peninsula would increase in the future. With CO 2 fertilization effects, Tao and Zhang 13 found that the relatively higher yield distributed in southwest of Shandong and south of Henan from 1961 to 1990 and could continue to increase in the future, while the relatively higher ET C distributed in southwest of Shandong and could continue to decrease in the future. The reason for these differences may exist in the GCMs, climate scenarios, downscaling methods, and crop models 13,14,18,32 . This study hypothesized that only the winter wheat cultivar of 'Bainong 207' would be sown and the planting boundaries would not change in Huang-Huai-Hai Plain. In fact, the wheat varieties and the planting boundaries will change, which in turn influence wheat yield and ET C . In the past three decades, new wheat cultivars were the dominated factor for yield increase in the North China Plain 35 . Water use efficiency increased substantially from 1.0-1.2 kg m −3 in early 1970s to 1.4-1.5 kg m −3 for recent cultivars 36 . In addition, agronomic practices and technological advancement could also affect ET C 5 . Therefore, there must be a great uncertainty in the prediction of water requirement of winter wheat in the future.
Some researches reported that the simulation accuracy was high under full irrigation and low under limited irrigation when the DSSAT model was used to simulate crop water requirement. Marek et al. 37 showed that the water requirement of maize simulated by DSSAT increased with the increase of deficit level and higher than the measured value. DeJonge 30 showed that the simulated water requirement of maize under full irrigation was slightly lower than the measured value, while it was much higher than the measured value under limited irrigation. Then he modified the model by integrating a dynamic K C , and found the simulation accuracy was improved for maize water requirement under both two irrigation conditions, especially under limited irrigation.

Conclusion
This study was to improve the DSSAT model to simulate ET C of winter wheat and understand the impacts of climate change on the potential yield and ET C of winter wheat in Huang-Huai-Hai Plain. The main conclusions drawn are as follows. The spatial distributions of potential yield in 2020s, 2050s, 2080s under RCP4.5 were similar to those under RCP8.5, they were yield increasing gradually from the northwest inland to the southeast coast. ET C was decreasing gradually from Shandong Peninsula to the surrounding area, and the minimum ETc was predicated in the southeast of Huang-Huai-Hai Plain in Jiangsu Province and the border between Henan and Anhui Provinces. The winter wheat potential yield, ET C and effective precipitation during the wheat growing period might increase from 2020s to 2080s under RCP4.5, while the IWR would decrease. Under RCP8.5, The effective precipitation during wheat season would decrease first from 2020s to 2050s, and then increase from 2050s to 2080s. The potential yield, ET C and IWR would increase first and then decrease.

Material and Methods
Data collection. Weather data. Daily meteorological data of 1961-2010 for 88 weather stations in Huang-Huai-Hai Plain were collected from the China Meteorological Data Sharing Service System (http://cdc.nmic.cn/ home.do) to calibrate the Hargreaves formula. The data included daily maximum and minimum temperatures, relative humidity, average wind speed and sunshine duration, and precipitation.
GCMs data. The Canadian Earth System Model-CanESM2 was used for future climate projection since it had been proved with a high accuracy in China 38 . The daily National Center for Environmental Prediction (NCEP) reanalysis dataset included twenty-six large-scale atmospheric variables used to calibrate and validate the SDSM model, including mean sea level pressure, near surface relative humidity, surface specific humidity, mean temperature at 2 m height. The same daily CanESM2 atmospheric variables under RCP4.5 and RCP8.5 at a resolution of 2.8125×2.8125° and for the period of 1961-2100, which were obtained from Canadian Climate Data and Scenarios (http://www.cccsn.ec.gc.ca), were used to generate future climate data.
DSSAT input data. Both field experiment and DSSAT model were used to evaluate the influences of climate change on the ET C of winter wheat in Huang-Huai-Hai Plain. The field experiment was conducted during the winter wheat growing seasons of 2015 and 2016 at the Experimental Station located in Xinxiang (35°18′ N, 113°54′ E, 73.2 m), Henan Province, China. Climate in this region is a temperate continental climate with a mean annual temperature of 14.1 °C and an average annual rainfall of 582 mm, the groundwater level is 8 m below the ground. A main cultivated winter wheat cultivar of 'Bainong 207' in northern Henan Province was sown in drills with a row spacing of 0.25m, depth of 0.05 m, and density of 400 seeds m −2 . There was no water and fertility stress, and insect and pest were well controlled during the whole growth periods. Observed ET C of each stage was calculated using a water balance method 39,40 (Eq. 1).
where ET C is the evapotranspiration of winter wheat (mm), M1 is the initial soil moisture (mm), P is the precipitation (mm), I is the irrigation (mm), CR is the water used by crop through capillary rise from groundwater (mm), and is negligible because the groundwater table is lower than 8 m below the ground surface 41 . M2 is the final soil moisture (mm), D is the deep drainage(mm), R is the runoff (mm). The division of each stage was based on the phenology of winter wheat, the date of irrigation and precipitation (the value of deep drainage and runoff we couldn't measure), and the measurement date of soil moisture. The divided stages were listed in Table 4. dates of the 88 stations were obtained from local investigation and the field management was set as the same as that of the field experiment.
The weather data of field experiment, including temperature, relative humidity, average wind speed, sunshine hour, and precipitation, was collected at 0.5 hours intervals from a weather station in Xinxiang. However, for the model simulation from 2011 to 2100 the weather data were projected by SDSM4.2, including maximum temperature, minimum temperature, precipitation, and solar radiation. In addition, the solar radiation from 1961 to 2005 was calculated by the Angstrom formula as follows (Eqs 2-6).
where R s is the solar short wave radiation (MJ m −2 d −1 ), R a is the solar radiation at the top of the atmosphere (MJ m −2 d −1 ), n is the actual sunshine duration (h), N is the maximum possible sunshine duration (h), a s , and b s are the regression constants, whose values change with the different meteorological conditions (e.g. humidity, dust) and solar declination (latitude, month).The suggested values of a s , and b s in Huang-Huai-Hai Plain are 0.152 and 0.556 in spring (from March to May), 0.115 and 0.588 in summer (from June to August), 0.301 and 0.311 in autumn (from September to November), and 0.172 and 0.536 in winter (from December to the following February) 42 . And dr is the reciprocal of the relative distance between the sun and the earth, ω s is the solar hour angle (rad), ψ is the geographic latitude (rad), δ is the solar declination (rad), J is Julian day of a year, D is the days of a year. The soil data for DSSAT simulation were partially obtained from the Harmonized World Soil Database V1.2 by the Food and Agriculture Organization of the United Nations (FAO) and the International Institute for Applied Systems Analysis (IIASA). The rest of soil data were got from the 1:1000000 Chinese Soil Map by the Nanjing Institute of Soil Science, Chinese Academy of Sciences. The database contained soil physical properties of two layers of topsoil (0-30 cm) and subsoil (30-100 cm), including sand, silt and clay fractions, soil bulk density, cation exchange capacity, organic carbon content and PH. The saturated soil water content, field capacity, and permanent wilting point of each soil layer were calculated by the following formula (Eq. 7-10) 43,44 .
where θ s , θ c , θ w are saturated soil water content, field capacity, and wilting coefficient, respectively, C is clay fraction, W is soil porosity, γ is soil bulk density (g cm −3 ), ρ is soil specific gravity, whose suggested value is 2.65 g cm −3 45 .  Methods. Statistical downscaling methods. At present, downscaling approaches include four main categories: regression methods, weather pattern approaches, stochastic weather generators, and limited-area regional climate models. Among these methods, statistical downscaling methods are the most widely used 46 . In this study, the statistical downscaling model SDSM4.2 was adopted for scale conversion. Following steps were involved to project future climatic variable. First, two sub-periods of 1961-1990 and 1991-2005 were chosen to calibrate and validate the model. Second, predictors were selected by seasonal correlation analysis, partial correlation analysis, and scatter plot analysis (α = 0.05). Third, the empirical statistical relationship between large-scale predictors and predictands (the maximum temperature, minimum temperature, solar radiation, and precipitation) was established to determine the parameters of multiple regression equation. And fourth, the coefficient of determination (R 2 ) and root mean square error (RMSE) were used to quantitatively assess the performance of SDSM. Finally, the girded data for the future climate provide by CanESM2 under RCP4.5 and RCP8.5 during three future periods of the 2020s (2011-2040), the 2050s (2041-2070), and the 2080s (2071-2100) were input into the models to generate the downscaled future daily climatic series for each station. where ET 0 is the daily reference crop evapotranspiration (mm d −1 ), EEQ is the equilibrium evapotranspiration (mm d −1 ), A LBEDO is the crop albedo, T max and T min are the daily maximum and minimum temperature, respectively (°C), S R is the daily total solar radiation (MJ m −2 d −1 ), α is the a coefficient of advectivity, which was set as 1.1 in the DSSAT V4.6 when temperature was between 5 and 35 °C and could be slightly smaller or larger than 1.1 (Eqs 13-14) when temperature was below 5 °C and above 35 °C.   In this study, the nonlinear regression analysis was used to calibrate the Hargreaves formula with the software of SPSS 31 . The main process was as follows. First, the daily maximum temperature, minimum temperature, and extraterrestrial radiation were chosen as independent variable and daily ET 0 calculated by the P-M formula were chosen as dependent variable for the 88 stations from 1961 to 2012. Next, the original values of K, n and T off of Hargreaves formula were set to 0.0023, 0.5 and 17.8 in SPSS. Third, new parameter values were obtained through iterative nonlinear regression analysis. Finally, correlation index and standard error were used to quantitatively assess of accuracy of model calibration.
Dynamic exponential decay function of crop coefficient (K c ). The DSSAT4.6 mode code employs the following formula for calculation of the crop efficient K c (Eq. 17). is still lower, for P-M method, the value of EORATIO is set to 1.1, which ensures K C varies daily between 1.0 and 1.1, and for PT method, K C is lacking. DeJonge et al. 30 suggested that the Kc was a dynamic exponential decay function of leaf area index (LAI) (Eq. 18).
C The verification of main parameters and modification plan of DSSAT. In this study, the parameter estimation and verification of the DSSAT-CERES-Wheat model was carried with the DSSAT-GLUE package. The main parameters involved were P1V, P1D, P5, G1, G2, G3, and PHINT. The estimation of parameter values needed to meet three conditions. First, phenological stages or simulated dates of flowering and maturity were consistent with the actual situation. Second, the accumulation speed of simulated biomass was consistent with the actual situation through adjusting the parameters related to leaf expansion and photosynthesis. Third, the simulated final yield of winter wheat was consistent with the actual yield. Due to the low accuracy of PT method for estimating ET C of winter wheat, the calibrated Hargreaves method and a dynamic exponential decay function of K c were adopted in the DSSAT4.6 model. Three DSSAT modification plans, including original PT model (Plan 1), the original PT model but with a dynamic K c function (Plan 2), and the calibrated Hargreaves model and a dynamic K c function (Plan 3), were used to simulate the growth of winter wheat under full irrigation in 2015 and 2016, respectively. Then a comparison between simulated ET C and observed ET C were made to select the modification plan with the highest accuracy. Finally, the potential yield, ET C and irrigation water requirements (IWR) of winter wheat in Huang-Huai-Hai Plain under RCP 4.5 and RCP8.5 were predicted by the selected one.
Spatial interpolation. The ordinary Kriging method was based on the spatial position between the measured and measuring points to make a linear unbiased optimal estimation for measuring points. Then the spatial distributions of factors in the study area were analyzed by generating a spatial distribution map for the key parameters. The detailed steps of ordinary Kriging in the ArcGIS v10.2 were as follows. First, input the basic geographic information (longitude, latitude, and elevation) and factor data into ArcGIS v10.2 to get the point layer for the 88 stations. Next, combined with the surface layer of Huang-Huai-Hai Plain boundary, the ordinary Kriging interpolation method in geostatistics was used to process the data. Finally, the values of the fitting parameters were checked by semi-variance or covariance functions, and the spatial distribution map of key parameters was generated.
i i i where P i is the i-th simulation value, Q i is the i-th observation value, P is the average of simulation value, Q is the average of observation value; and n is the number of samples.