agronomy The Adaptability of APSIM-Wheat Model in the Middle and Lower Reaches of the Yangtze River Plain of China: A Case Study of Winter Wheat in Hubei Province

: The middle and lower reaches of the Yangtze River (MLYR) plain represent the second-largest wheat producing area in China; the winter wheat-rice system is one of the main planting systems in this region. The use of the agricultural production system simulator (APSIM)-wheat model to simulate wheat production potential and evaluate the impact of future climate change on wheat production in this region is of great importance. In this study, the adaptability of the APSIM-wheat model in the MLYR was evaluated based on observational data collected in ﬁeld experiments and daily meteorological data from experimental stations in Wuhan, Jingmen, and Xiangyang in Hubei province. The results showed signiﬁcant positive relationships between model-predicted wheat growth duration from sowing to anthesis and maturity and the observed values, with coe ﬃ cients of determination (R 2 ) in ranges of 0.90–0.97 and 0.93–0.96, respectively. The normalized root-mean-square error (NRMSE) of the simulated growth durations and measured values were lower than 1.6%, and the reﬁned index of agreement (d r -values) was in the range of 0.74–0.87. The percent mean absolute relative error (PMARE) was cited here as a new index, with a value below 1.4%, indicating that the model’s rating was excellent. The model’s performance in terms of grain yield and above-ground biomass simulation was also acceptable, although it was not as good as the growth periods simulation. The R 2 value was higher than 0.75 and 0.72 for the simulation of grain yield and biomass, respectively. The indices NRMSE and PMARE were lower than 19.8% and 19.9%, and the d r -value was higher than 0.71. According to our results, APSIM-wheat was an e ﬀ ective and accurate model for simulating the phenology and yield production processes of wheat in the MLYR, and the results also provided a theoretical basis and technical support for further research on the yield potential of wheat-rice rotation planting systems with clariﬁcation of the key factors limiting the yield gap in this region.


Introduction
In China, more than 50% of the population is mainly fed by wheat, and with less and less arable land being available, further improvement in wheat production per area is crucial to ensure China's

Experimental Site
Hubei province (29 • 01 -33 • 6 N, 108 • 21 -116 • 07 E) is located in the central part of China, and its environment is representative of the typical conditions in the MLYR plain. The area exhibits a subtropical monsoon humid climate with annual average sunshine hours of 1754-1910 h, an annual average accumulated temperature (≥0 • C) of 5266-8372 • C d, abundant precipitation in the wheat-growing season, and annual average precipitation of more than 1000 mm. These are superior agricultural production conditions and are conducive to achieving a good crop yield. In this paper, three winter wheat production regions were selected (Figure 1), namely, Wuhan (WH), Jingmen (JM), and Xiangyang (XY), where the winter wheat growth areas account for over 70% of the total wheat planting area in Hubei province. Figure 1. Locations of the experimental sites (black solid dots) that were used to calibrate and validate the agricultural production system simulator (APSIM)-wheat model. Note: The experimental sites were Guyi town in Xiangzhou county, Zhangjiaji town in Xiangzhou county, Zaoyang county, and Yicheng county in Xiangyang; Zhongxiang county, Qujialing county, and Shayang county in Jingmen; Wuchang district and Hongshan district in Wuhan.

Data Source
In terms of weather data, records of daily maximum and minimum temperatures, sunshine hours (or daily solar radiation), and precipitation were available from 1999 to 2016 at each experimental weather station and the National Meteorological Information Center (CMDC, http://data.cma.cn/en). As actual solar radiation measurements were not available for most of the stations, the daily solar radiation was estimated from the sunshine hours using the Angstrom formula [20].
For experimental data, information on winter wheat phenology aspects (such as sowing, emergence, jointing, anthesis, maturity dates, cultivar type, yield, and biomass) and management practices (such as planting density and fertilizer) was obtained from each experimental site (previous field experiment or local agricultural technology extension center) and from the literature (listed in Table 1). Management details of each experimental site can be found in Supplementary Materials (Tables S1 and S2). Crop management practices at each experiment site were generally the same as, or better than, the local traditional practices [21]. The wheat varieties planted at each site are also listed in Table 1. These were widely used in each province during the planting years.
Regarding soil data, the soil bulk density (BD), saturated volumetric water content (SAT), drained upper limit (DUL), 15-bar lower limit (LL15), soil organic carbon (SOC), and pH value (pH) in different soil layers were obtained from the China Soil Scientific Database (CSSD, http://www.soil.csdb.cn), the dissertation of Sun [22] and our previous measurements. Table 1. Wheat varieties and sources of experimental data for agricultural production system simulator (APSIM)-wheat model calibration and validation in Hubei province.

APSIM-Wheat Model
The APSIM-wheat model used in this study includes the crop-wheat, soil, surface organic matter (OM), and manager modules. The model operates on a "daily" time step and simulates daily crop development, biomass production, soil moisture, and nitrogen dynamics as affected by the climate and management measures.
Crop ontogeny aspects, such as daily biomass and yield production, are simulated by the relationships among crop growth, temperature, and photoperiod [33]. Crop biomass accumulation depends on the solar radiation interception (RI) and radiation use efficiency (RUE), taking full account of the effects of water and nitrogen constraints. RI and RUE related to the leaf area in this model, and the leaf area is calculated by the increase in the leaf dry weight and the maximum specific leaf area (SLAmax), which is related to the leaf area index (LAI), so a function (SLAmax = hSLA (LAI)) is defined by the parameters x_lai and y_sla_max in the APSIM-wheat model [12]. Crop development is primarily based on thermal time, whereas leaf and stem growth rates are calculated depending on phenological stages. The soil water module includes the lower limit (LL15), the drained upper limit (DUL), and the saturated (SAT) volumetric water contents of a sequence of soil layers. Further detailed descriptions of the APSIM-wheat model structure and processes were described by Keating et al. [12] and are available at the following website: http://www.apsim.info.

Calibration and Validation of the APSIM-Wheat Model
The data for the winter wheat growth period, leaf area index, above-ground biomass at anthesis and maturity, grain yield, yield components, winter wheat management at all stages, soil water, and nutrients content at a soil depth of 0-240 cm with 15 cm per layer were collected from the Wuhan, Jingmen, and Xiangyang experiment regions from 1999 to 2016. The database was divided into two parts (Table 1): one was used to calibrate the genetic parameters of the four wheat varieties (Zhengmai9023, Emai596, Emai18, and Emai170), and the other one was used to validate the utility and adaptability of the calibrated model.
In this study, the crop variety parameters for calibration, such as sensitivity to vernalization, sensitivity to the photoperiod, the thermal time required from grain filling to maturity, the number of grains per gram of stem, the potential grain filling rate, and the maximum specific leaf area, are shown in Table 2. A trial-and-error method was used to identify the optimal parameters of the three varieties. The growth periods from sowing to anthesis and sowing to maturity were mainly calibrated by the sensitivity to vernalization, the sensitivity to the photoperiod, and the thermal time required from grain filling to maturity. The grain yield and above-ground biomass were mainly calibrated via the number of grains per gram of stem, the potential grain filling rate. A maximum specific leaf area was used to adjust the leaf area index (LAI) during the growing season. Note: † . ∆LAI means delta leaf area index.

Evaluation of the Model's Performance
The coefficient of determination (R 2 ), root mean square error (RMSE), normalized root mean square error (NRMSE), Willmott's refined index of agreement (d r -value) [34], and percent mean absolute relative error (PMARE) [35] were used to evaluate the utility and adaptability of the calibrated model. R 2 can reflect the consistency between simulated and observed values-the closer to 1, the better the performance is. The RMSE and NRMSE reflect the relative error and absolute error between simulated and observed values-the smaller the value, the better the performance is. The d r -value (in the range of −1 to 1) indicates the sum of the magnitudes of the differences between the model-predicted and observed deviations of the observed mean relative to the sum of the magnitudes of the perfect-model and the observed deviations of the observed mean [34]. The PMARE directly indicates the strengths or weaknesses of the simulation and thus helps to decide whether to accept or reject the model. The suggested performance rating for model evaluation based on the PMARE is listed in Table 3 [35]. Table 3. Suggested performance rating for model evaluation based on the percent mean absolute relative error (PMARE) [35].

PMARE Value (%)
Model Rating The equations are listed below: where P i and O i are the predicted (or simulated) and the observed value, respectively; O is the mean of the observed values; n is the number of observed values.

Statistical Analysis of the Data
The weather (including daily maximum and minimum temperatures, sunshine hours (or daily solar radiation), and precipitation) and crop management (such as sowing date, irrigation, nitrogen application, and so on) data were compiled in Excel (Microsoft Office Professional Plus 2010, Microsoft Corporation, WA, USA) and then used to calibrate the model and simulate the phenology, grain yield, and above-ground biomass of each variety at the three experimental sites. Regression analyses between simulated and measured values were performed using, t 12.5 software (Systat Software, Inc, San Jose, CA, USA).

Genetic Parameters of Each Variety for the APSIM-Wheat Model
A "trial-and-error" method was used to determine the genetic parameters of Zhengmai9023, Emai596, Emai18, and Emai170 for the APSIM-wheat model ( Table 4). The sensitivity to vernalization and the photoperiod ranged from 2.4 to 3.5 and 2.5 to 3.6, respectively, with the values from Zhengmai 9023 being the highest. The thermal time required from grain-filling to maturity of Zhengmai 9023 was 650 • C d, which was higher than that of the other varieties, and the value from Emai596 was the lowest, accumulating only 520 • C d. The number of grains per spike was over 28 for all varieties, and the potential grain filling rate reached to 0.0031 g grain −1 d −1 . The maximum specific leaf area ranged from 21,000 to 25,000 mm 2 g −1 when ∆LAI (x_lai = 0), and from 18,000 to 23,000 mm 2 g −1 when x_lai = 5. Table 4. Genetic parameters, sensitivity to vernalization (STV), sensitivity to photoperiod (STP), thermal time required from grain-filling to maturity (TTRGM), numbers of grain per stem (NGS), potential grain filling rate (PGFR), and maximum specific leaf area for ∆LAI (MSLA) for the model simulation of the four wheat varieties in Hubei province.

Performance of the Calibrated APSIM-Wheat Model
Wheat growth duration, grain yield, and above-ground biomass data were used to validate the performance of the calibrated APSIM-wheat model. Relative data came from the "validation part" shown in Table 1. The data from the "calibration part" were not included because they were used to calibrate the model, but they can be found in the Supplementary Materials ( Figure S1 and Table S3).

Wheat Growth Duration
The simulated wheat growth duration in Wuhan, Jingmen, and Xiangyang was significantly positively correlated with the measured values ( Figure 2 and Table 5), with the coefficient of determination (R 2 ) for the duration from sowing to anthesis being over 0.90. The NRMSE values between the simulated and observed values of the experimental sites in Wuhan, Jingmen, and Xiangyang were low at 1.2%, 1.3%, and 1.6%, respectively, and the d r -values were high at 0.81, 0.78, and 0.79, respectively. The R 2 values between the simulated and observed durations from sowing to maturity ranged from 0.93 to 0.96, with the NRMSE values being 1.1%, 1.0%, and 1.3%, and d r -values being 0.87, 0.76, and 0.74, in Wuhan, Jingmen, and Xiangyang regions, respectively. The PMARE values ranged from 0.8 to 1.4, in the range of an excellent model rating. According to the comparison between the simulated results and the measured values, the mean difference in wheat growth duration at the three test regions was within three days (Table 6). This showed that the APSIM-wheat model provided an excellent simulation performance during the determination of the growth duration of winter wheat in the middle and lower reaches of the Yangtze River plain of China.

Grain Yield and Above-Ground Biomass
The crop grain yield and above-ground biomass were the key indices used to evaluate the utility of the suggested model. Significant positive relationships between the simulated and observed values of wheat grain yield and above-ground biomass at three experimental regions in Wuhan, Jingmen, and Xiangyang were observed with the R 2 values of 0.75-0.78 for the grain yield 0.72-0.87 for the above-ground biomass (Figure 3 and Table 7). The NRMSE values for the simulated and observed values of grain yield were 13.0%, 18.1%, and 17.1%, and the d r -values were 0.73, 0.78, and 0.71 in the above experimental regions. The NRMSE values for the simulated biomass and the observed values were 13.5%, 17.7%, and 19.8%, and the d r -values were 0.80, 0.71, and 0.72 in the above experimental regions. Thus, the NRMSE values for the wheat grain yield and above-ground biomass in the three trial regions were all less than 20%, and the d r -values were all above 0.71, and the mean differences between the simulated and observed wheat grain yield and above-ground biomass were lower than 556 and 1323 kg ha −1 , respectively (Table 6). Meanwhile, the PMARE values ranged from 10.6% to 19.9% in most of the model equations performed, which were in the range of either good or fair (Table 3). This indicated that the APSIM-wheat model could simulate the grain yield and biomass of winter wheat in the middle and lower reaches of the Yangtze River plain well and could effectively and accurately simulate the wheat yield formation process in this region.  Table 7. The fitted equations for the linear relationships between the model-simulated and observed wheat grain yield and above-ground biomass and reliability measures: coefficient of determination (R 2 ), normalized root mean square error (NRMSE), refined index of agreement (d r -values), and percent mean absolute relative error (PMARE).

Discussion
The APSIM model simulates the multi-parameter interactions of soil, crops, meteorology, etc., but it takes soil as the center and fully considers the impacts of intercropping and rotation on the soil. It has been well applied in different planting systems in Australia, America, The Netherlands, New Zealand, Germany, and other parts of the world [13][14][15][16]. Seyoum et al. used the APSIM model to examine G × E × M interactions for maize improvement in Ethiopia and found that the model accurately predicted plant-available soil water, days-to-flowering, days to maturity, the leaf area index, biomass, and the yield of maize [36]. Makowski et al. concluded that the APSIM-wheat model could predict phenological stages, like anthesis and maturity, at values close to the observed values [37]. In China, the APSIM-wheat model has been used in the Northeast, Northern, Northwest, and Southwest regions of China to simulate crop growth processes and evaluate the effect of climate change on crop production. This should indicate that the APSIM-wheat model could be used to reproduce the observed crop growth, yield, potential yield, and water use in the study areas [38,39]. The MLYR plain is located in the middle of China and is the second-largest wheat producing area. Wheat production in this region is risky because of the high climatic variability in recent years. Therefore, effective and accurate models are urgently needed to study the wheat growth, development, and yield in response to climate variability in this region.
The APSIM-wheat model was used for the first time and evaluated winter wheat in wheat-rice continuous cropping systems in the MLYR plain in this study. Accurate phenology is the priority when calibrating crop models [40]. The phenology of wheat has a strong influence on the development and grain yield of the crop [41]. Simulation outcomes from the APSIM-wheat model in our study showed that this model could be used as a suitable tool for the selection of appropriate cultivars and to investigate the effect of climate variability on wheat growth and yield [11]. Li et al. found a good agreement between values simulated by the APSIM-wheat model and the observed values for the flowering and maturity dates, yield, and biomass in the North plain of China, with R 2 values in the range of 0.72-0.86, the d-values of 0.90-1.00, and RMSE values of 1.10-3.95 d, respectively, for the observed growth periods. The simulated yields were close to the observed yields, with R 2 > 0.79, d-values > 0.90, and NRMSE < 13.4%, respectively [39]. We used the same genetic parameters of Zhengmai9023, which were obtained by Li et al. [39] in the North plain of China to validate the APSIM-wheat performance in Xiangyang (to represent MLYR). The poor performance of APSIM-wheat was observed for the simulation of the growth to maturity duration (with 187 d and 216 d for the simulated and observed values, respectively) and grain yield (with 2965 and 7936 for the simulated and observed values, respectively). This indicated that the genetic parameters of Zhengmai9023 from the North plain China could not be used in the model to simulate yield production in the MLYR. Therefore, it was important to calibrate and validate the APSIM-wheat model before it was adopted in a new environment or region. In this study, the calibrated APSIM-wheat model was shown to predict the phenological stages well in the MLYR plain of China. For instance, the predicted values for the sowing date to anthesis and sowing to maturity were close to the values of observed data. The R 2 values were all over 0.90, and the NRMSE values were all under 1.6% (Figure 2 and Table 5). The results of the simulations conducted with this model showed that the yield and biomass values were also close to the observed values ( Figure 3, Tables 6 and 7), with R 2 > 0.72, d r -values > 0.71, and NRMSE < 19.8%, respectively. These results showed that process-based models had good potential to simulate crop biomass and yield production in the MLYR [42].
The use of logical, consistent, and generally accepted indicators is important when discussing the evaluation of the model procedure. The indicators should appropriately quantify the objectives of the model evaluation and indicate its utility. Fox recommended that the mean error, mean absolute error, variance of the distribution of difference, and root mean square error (R) or its square (R 2 ) should be calculated and reported when evaluating the performance of models [43]. Willmott commented that the correlation between model-predicted and observed data, commonly described by Pearson's product-moment correlation coefficient, that is R and R 2 , is an insufficient and often misleading measure of accuracy [34]. Ali et al. also reported that the difference-based (R, R 2 , RMSE, NRMSE), efficiency-based (Nash and Sutcliffe coefficient, model efficiency of Loague and Green, Legates and McCabe's index) measures were found ambiguous, inconsistent, and not logical in many cases [35]. Thus, in addition to R 2 and NRMSE, two new indices (d r -value and PMARE value) were cited and calculated in this study. The d r -value is intended to be a descriptive measure, and it is both a relative and bounded measure that can be widely applied in order to make cross-comparisons between models [34]. It was shown that the d r -values of the model-predicted wheat growth periods versus the measured values ranged from 0.74 to 0.87, with a mean value of 0.79, indicating the good performance of the model in terms of predicting the wheat growth period. Ali et al. also observed that PMARE always followed logical behavior and no ambiguous result [34]. The PMARE value could be used to evaluate model performance. For the growth periods of anthesis and maturity, the PMARE values indicated an excellent rating in this study (Tables 3 and 4). In addition, the model's performance for grain yield (PMARE < 14.6%) and biomass (PMARE < 19.9%) produced values that were acceptable (Tables 3 and 7) [35], although they were not as good as the growth period simulation.
Further analysis indicated significant positive relationships between the indicators PMARE and NRMSE, with the R' 2 of 0.96 for the simulation values. Significant negative relationships were observed between PMARE and R 2 and between NRMSE and R 2 , with the R' 2 of 0.67 and 0.76, respectively. The index d r -value had no clear consistent correlation with the three other parameters, and this result was consistent with Willmott et al. [34] (Figure S2).
Quantitative predicting of the phenological period is of great significance for accurately simulating crop growth and productivity formation [44]. There are 11 phases in APSIM-wheat model, and the timing of phases after germination is determined by the accumulation of thermal time, which is adjusted by vernalization, photoperiod, and other factors, such as nitrogen, water. The length of each phase is determined by the thermal time target, and the next phenological stage occurs when an adjusted thermal time reaches the target thermal time. In the APSIM-wheat model (7.5 R3008 version), soil, water, nitrogen, and phosphorus stresses have no effect on phenological development, so the adjusted thermal time is determined by daily thermal time, vernalization, and photoperiod factors [45]. Thus, we decided to calibrate vernalization and photoperiod to predict these phases in this study, the same as some previous studies [39,46]. As mentioned above, the next phase will begin once an adjusted thermal time reaches the target thermal time. The target thermal time is a site-cultivar-special parameter, and we did not calibrate tt_end_of_juvenile and tt_floral_initiation in this study and only calibrated vernalization and photoperiod-related parameters for days so that anthesis would increase the chances of a less robust calibration and hence a lot of uncertainty in the model implementation, which would be explored in our next step.
Leaves are important organs for intercepting solar radiation and then producing dry matter for the whole plant. Simulating leaf area index (LAI) is one of the important output parameters for a crop growth simulation model [47,48]. In the APSIM-wheat model, the daily increase in leaf area is the minimum between stressed LAI and the carbon-limited LAI. The stressed LAI is determined by leaf number, leaf area, plant population, and environments, such as nitrogen, soil, water, and phosphorus stress factors. Leaf area related to carbon production is calculated by the increase in dry weight of leaf and the maximum special leaf area (sla-max); the sla-max is the maximum leaf area that can be expanded per gram of biomass, which regulates leaf area index and further regulates radiation interception, generating biomass with RUE [45]. In this study, we calibrated the sla-max to adjust LAI for wheat. The validation results showed that the simulated LAIs were close to observed values, the R 2 was in the range of 0.85-0.98, NRMSE being from 6.9% to 12.5%. The slopes of the trend in LAI-observed over LAI-simulated for jointing, anthesis, and the whole data were 1.10, 0.93, and 0.96, respectively (Table S4 and Figure S3). This result indicated it was justifiable by using the values of sla-max in Table 4 to calibrate LAI. We also simulated LAIs at jointing and anthesis with the default values of sla-max in the APSIM-wheat model (Table S4 and Figure S3), and the results showed that the LAIs at jointing and anthesis were overestimated (slopes of a trend in LAI-observed over LAI-simulated were in the range of 0.37-0.48). In addition, the NRMSE was in the range of 28.6-36.5%, even though the p-values for the fitted lines were below 0.01, indicating a poor performance of APSIM-wheat on LAI with the default values of sla-max. Anyway, considering the calculation equation for LAI, it would be better to collect more data, including leaf number and leaf area, to improve the performance of the APSIM-wheat further in the next step.
A major source of uncertainty in our analysis may derive from the trial sites and cultivars chosen to calibrate the model [38]. In this study, only Hubei province was selected to test the performance of the APSIM-wheat model in the MLYR plain of China. Although this province has the typical subtropical climate conditions of this region, it still has its limitations in representing the varied topography and cropping systems of the region. Meanwhile, four representative cultivars, which are the main cultivars planted in Hubei province, were chosen as the experimental data to calibrate the model. More than 340 cultivars were released and introduced into the national seed market in China from 1984 to 2010 [39], and approximately 40% were planted in the MLYR plain. Another source of uncertainty may derive from the data from wheat trials coming from the literature and local agricultural technical departments. The data from the field experiment generally ignored the effects of diseases, insect pests, and weeds. This might have influenced the performance of the APSIM-wheat simulations. Therefore, in future research, we would expand the selection range of trial sites and wheat cultivars to further optimize the combination of model parameters and to improve the accuracy and reliability of the model. Then, it would be possible to use the validated model to simulate the yield potential in this region, and the researchers might be able to clarify the potential yield, yield gap, and key limiting factors across the MLYR and develop or optimize wheat production management to further improve grain yield. Furthermore, this model could be used to explore the wheat yield performance under conditions of future climate change, which would help breeders to select and breed varieties that could adjust to climate change. It would also help policymakers to develop methods to optimize the distribution of crops and varieties to ensure national food security.

Conclusions
The APSIM-wheat model was used to predict the wheat potential grain yield, yield gap, and yield production response to climate variability in the middle and lower reaches of the Yangtze River plain of China. The utility and adaptability of the model were evaluated and discussed in this study. Significant positive relationships were observed between the simulated and observed wheat growth durations, grain yield, and above-ground biomass. Meanwhile, some logical, consistent, and generally accepted indicators-the coefficients of determination, normalized root-mean-square error, refined index of agreement, and percent mean absolute relative error-were used to evaluate the performance of this model, and the results indicated that the APSIM-wheat model had an excellent performance for determining wheat growth periods and an acceptable performance for determining grain yield and above-ground biomass, although its performance was not as good as the growth period simulation. According to our results, using the APSIM-wheat model was an effective and accurate way to simulate the phenology and yield production process of wheat in the MLYR. The results also provided a theoretical basis and technical support for further research on the yield potential of a wheat-rice planting system and to clarify the key limiting factors related to the yield gap in this region.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/7/981/s1, Figure S1: Relationships between simulated and observed growth duration from sowing to anthesis (A) and to maturity (B), grain yield (C), and above-ground biomass (D) of winter wheat with calibration dataset in Hubei province, Figure S2: Relationships between indices percent mean absolute relative error (PMARE), determination coefficients (R 2 ), normalized root-mean-square error (NRMSE), and index of agreement (dr-value) for all simulations, Figure S3: Relationships between simulated and observed LAI at jointing and anthesis with sla-max values from Table 4 (a) and default sla-max in the APSIM-wheat model (b), Table S1: Details of nitrogen application rate and plant density of each experimental site for APSIM-wheat model calibration in Hubei province, Table S2: Details of nitrogen application rate and plant density of each experimental site for APSIM-wheat model validation in Hubei province, Table S3: The fitted equation for the linear relationships between model-simulated and observed growth duration, grain yield, and biomass with calibration dataset and reliability measures: coefficient of determination (R 2 ), normalized root mean square error (NRMSE), refined index of agreement (d r ), and percent mean absolute relative error (PMARE), Table S4