CROPGRO-Soybean Model Calibration and Assessment of Soybean Yield Responses to Climate Change

Process-based crop simulation models are useful for simulating the impacts of climate change on crop yields. Currently, estimation of spatially calibrated soil parameters for crop models can be challenging, as it requires the availability of long-term and detailed input data from several sentinel sites. The use of aggregated regional data for model calibrations has been proposed but not been employed in regional climate change studies. The study: 1) employed the use of county-level data to estimate spatial soil parameters for the calibration of CROPGRO-Soybean model and 2) used the calibrated model, assimilated with future climate data, in assessing the impacts of climate change on soybean yields. The CROPGRO-Soybean model was calibrated using major agricultural soil types, crop yield and current climate data at county level, for selected counties in Alabama for the period 1981-2010. The calibrated model simulations were acceptable with performance indicators showing Root Mean Square Error percent of between 27 43 and Index of Agreement ranging from 0.51 to 0.76. Projected soybean yield decreased by an average of 29% and 23% in 2045, and 19% and 43% in 2075, under Representative Concentration Pathways 4.5 and 8.5, respectively. Results showed that late-maturing soybean cultivars were most resilient to heat, while late-maturing cultivators needed optimized irrigation to maintain appropriate soil moisture to sustain soybean yields. The CROPGRO-Soybean phenological and yield simulations suggested that the negative effects of increasing temperatures could be counterbalanced by increasing rainfall, optimized irrigation, and cultivating latematuring soybean cultivars.


Introduction
Global climate change and variability is expected to affect crop phenological processes and yield. The need for the assessment of the extent of such potential impacts have resulted in the use of several process-based crop models to simulate the unique interactions between crops, their environment (weather and soils), management practices and performance under projected future climate conditions (Ahmed & Hassan, 2011). One commonly used process-based crop model is the Crop Growth (CROPGRO) (Jones et al., 2003). It is a grain legume model based on the SOYGRO, BEANGRO and PNUTGRO models (Tsuji et al., 1998). Numerous calibration and validation studies on CROPGRO-soybean have established the model's ability to simulate crop development and observed seed yield under different climate conditions and regions (Batchelor et al., 1993;Tsuji et al., 1998;Lal et al., 1999;Southworth et al., 2002;Wang et al., 2003;Mera et al., 2006;Res et al., 2007).
The successful use of these crop simulation models depends on the availability and use of long term experimental data for accurate calibration (Lobell & Ortiz-Monasterio, 2006;Cabrera et al., 2007;Liu et al., 2011). However, the calibration of these models in regional climate change studies has numerous challenges. First, there is difficulty in obtaining long-term (30 or so years) data of detailed and complete regional crop growth and development information, crop yields and crop management records for use in conventional model calibration (Jiang et al., 2014). Second, soil survey reports fall short in their descriptions of soil characteristics, especially and with respect to details on physical and chemical parameters needed to approximate soil components (water, air, nutrients) and soil structure which affects root distribution (Mavromatis et al., 2001;. One way to overcome these challenges is using detailed data (e.g. time-series data on phenology, growth, soil nitrogen status, etc. and/or end-of-season yield and yield-component data) obtained from sentinel sites such as agricultural experiment stations and representative farms. Through this approach, several studies have successfully assessed the impact of climate change at the regional scale. Alexandrov & Hoogenbom (2000) calibrated and validated their CERES models for Maize and Wheat using detailed experimental data from a variety of trial sites across Bulgaria for the period 1980-1993. In their assessment of climate change and climate variability impacts on corn yields in the midwestern US, Southworth et al. (2002) used detailed experimental data for 1975-1990 obtained from two selected representative farms in Illinois to calibrate and validate CERES-Maize models (Southworth et al., 2002).  used experimental data for 1961-1990 from agricultural research stations to calibrate CERES, CROPGRO and SUBSTOR models in order to study climate change impacts on crops in southern Québec, Canada.
A second approach is estimating spatially variable soil parameters in a systematic, stepwise, and subjective approach for the calibration of crop models.

Study Area
The study sites include four counties within the State of Alabama in the US, which are among the highest soybean producers in the state, and whose historical yield statistics were readily available to the public. Counties were also selected to represent a spatial distribution or trends across the state to better incorporate the different weather patterns between the southern, central, and northern parts of the state. The research Counties selected are Limestone

DSSAT Crop Simulation Model: CROPGRO-Soybean Model
The Decision Support System for Agrotechnology Transfer (DSSAT), v4.6, is a software suit of over 28 crop simulation models (Hoogenboom et al. 2015). The software application is supported by database management programs for soil, weather, and crop management and experimental data, and by utilities and application programs in its simulations of growth, development and yield as functions of the soil-plant-atmosphere dynamics (Hoogenboom et al. 2015). In this study, CROPGRO-Soybean model was used to simulate crop growth. CROPGRO-Soybean model simulates growth, development and yield of a crop using prescribed or simulated management together with changes in soil water, carbon and nitrogen that take place during a cropping system.

Weather Data
Daily weather data (minimum and maximum air temperature, rainfall, solar radiation) for each site were obtained from the US Department of Agriculture-Agriculture Research Service (USDA-ARS) website (http://ars.usda.gov/Research/docs.htm?docid=19422) and National Aeronautics and Space Administration-Prediction of Worldwide Energy Resource (NASA-POWER) websites (http://power.larc.nasa.gov/cgi-bin/cgiwrap/solar/agro.cgi?email=agroclim@larc. nasa.gov).  et al., 1973;Huang et al., 2015). The need for this procedure stems from the fact that technology influences positive yield trends for almost all crops in the US, including soybeans (Easterling et al., 1996;Huang et al., 2015;USDA-NASS, 2015). In this study, yields were adjusted using the linear trend shown below:

Crop Yield Data
where: i = year; j = county; ij Y ′ = is the adjusted crop (soybean) yield in the county (j) in year (i); ij Y = is the crop (soybean) yield for a county (j) in year (i); A = is the year to which the yields were adjusted/de-trended (i.e., 2004).

Crop Management Data
"Planting windows" were estimated from planting dates obtained from US Department of Agriculture National Agricultural Statistics Service (USDA-NASS) website (Table 1). The model was set to simulate planting when the lower limit of percentage of soil water reaches or exceeds 70 and the soil warms up sufficiently. The nitrogen application rate and depths as well as row and plant spacing, and depth used were based on recommendations by the Alabama Agricultural Experiment Stations (AAES). The study assumed irrigation and pesticide application were absent from the study sites and therefore were not considered.

Cultivars and Genetic Coefficients
The crop cultivars were characterized by a specific set of genetic coefficients the model uses to simulate and predict the daily growth and development of the plants in response to weather, soil conditions and management practices. The study employed the model's validated generic soybean coefficients for the different soybean Maturity Groups (MG) (Tsuji et al., 1994). The MG 5 was used for Limestone and De Kalb counties while MG 6 and MG 7 were used for Dallas and Baldwin counties, respectively (Table 2).

Soil Data
The majority soil series used for agricultural purposes were identified from USDA soil surveys and published works and used as representative soils for the study sites (Table 3).

Future Climate Scenarios
MarkSim DSSAT weather generator (http://gisweb.ciat.cgiar.org/MarkSimGCM/) was used to downscale future climate data from two General Circulation Models (GCM): MIROC 5 (Watanabe et al., 2010) and IPSL-CM5A-MR (Dufresne et al., 2013). At each study site, climate scenarios were generated for two future periods: "2045" (2030-2060), and "2075" (2060-2090). Future climatic conditions according to the Representative Concentration Pathway (RCP) 4.5 and 8.5  respectively, for medium and high radiative forcing conditions, were used for the study. The two study periods allowed for the assessment of climate change impact in the near term ("2045") and in the long term, ("2075"

Carbon Dioxide Fertilization Effect-DSSAT Simulations
Another aspect of climate change which impacts crop yields is carbon dioxide. Elevated carbon dioxide concentrations affects yields in two ways: promotion of photosynthesis and secondly through improved water use efficiency (Justino et al. 2013). To simulate the fertilization effect of CO 2 on crop physiology, carbon dioxide concentration levels were introduced directly in the environmental modifications of DSSAT (Justino et al., 2013).

Soil Parameterization
This study employed a hindcast method to calibrate selected soil parameters using stepwise procedure adopted from Mavromatis et al. (2001). The soil parameters which defined water availability, water holding capacity, fertility and root growth were selected, and their initial values were set for at field-based measured values for Alabama region according to Ratliff et al. (1983). The purpose of this was to have a good starting point to help optimize the calibration the soil parameters. The best values of the parameters were chosen by simulating 30 years historical yields in CROPGRO-Soybean model and using scatter plots of observed, which had been adjusted to remove technology effects and simulated yields to evaluate the performance of the models. The set of parameters that gave the highest d-statistic and the lowest RMSE was then chosen. To further assess the accuracy of the CROPGRO-Soybean model calibration and statistics, the study employed scatter plots to evaluate model simulations and measured data (Liu et al., 2011). The scatter plots of observed versus simulated yields allows the modeler to examine how the manual shifts of the different soil parameters affect yield simulations. The initial soil profile conditions used before the calibration process started were based on default model variables values and values reported by Ratliff et al. (1983), which were set to start on January 3 of each year, a day before the start of simulations. For each soil within the study sites, the Soil Drained Upper Limit (SDUL), Soil Lower Limit (SLL) of extraction, Upper Limit Saturated (SAT), Soil Fertility Factor (SLPF) and the Soil Root Growth Factor (SRGF) were parameterized. For each of the soils, the delta (difference between drained upper limit (DUL) and lower limit (LL)) for each of the layers in the profile was initially set to a similar starting point: 0.100 cm 3 /cm 3 . The delta was shifted by adding or subtracting 0.005 cm 3 /cm 3 to either the LL or DUL while checking the lower end of the scatter plots to see how well the simulated yields mimicked the observed yields in the water limited years of the baseline. Once shifts in the delta stopped improving model performance, attention was turned to the SAT values. The SAT values were then successively shifted by ±1 cm 3 /cm 3 while checking how simulated yield mimicked the observed yields in the scatter plots in the water-limited years. This study ensured that the SAT values were always greater than the DUL values as pointed out by Romero et al. (2012). It is worth noting that the initial SAT value of the top layer was adopted in all the successive layers before systematic shifts were made. SLPF (scale of 0 -1) for each of the soils was initially set to 1.00 to give a similar starting point and then shifted by ±0.01 while checking the higher end of the 1:1 scatter plots to see how well the simulated yields mimicked the observed yields in the good (high) yielding years in the baseline. After the optimization was complete, the SRGF was set to 1 in the layers whose center was ≤30 cm from the top of the surface layer. The initial starting point SRGF for the remaining layers beneath was estimated using the relationship established by Gijsman et al. (2007), and given as 1 × exponential (−0.02 Layer-Centre) i.e. depth from the top of the soil surface to the center of the layer of interest.

Climate Change Impact Assessment
To assess the response of soybeans to future climate scenarios, the study utilized the seasonal analysis program of DSSAT v4.6 whereby seasonal experiments were set up, model simulations were run, and finally biophysical analysis of model results was done. The crop management variables, genetic coefficients and optimized soil profiles from current climate simulations were taken as input data during the assessment. In addition to the current weather data, future climate scenarios obtained from MarkSim were included in the analysis. Comparisons were made between mean yields of the baseline and future climate scenarios, with and without carbon fertilization. To assess the effect of carbon fertilization on soybean, carbon dioxide concentration was changed in the environmental modifications section of DSSAT v4.6 in each of the future climatologies treatments, as shown in Table 4.

Adjusted Yields
The de-trending/adjustments of soybean yields resulted in trend line gradients close to zero for all counties (Figure 3). The remaining residuals in the graphs indicate the inter-annual variations in yields due to weather (Easterling et al., 1996;Osborne & Wheeler, 2013;Huang et al., 2015). The adjustment resulted in soybean yields decreases to attain a period of approximately constant technology.

Model Calibration and Estimated Soil Parameters
The Soil Fertility Factor (SLPF) values ranged from 0.572 to 0.758 (Table 5), which is within the range of values supported in the literature Mavromatis et al., 2001). Delta (DUL-LL) ranged from 0.085 cm 3 /cm 3 to 0.220 cm 3 /cm 3 for the soil profile layers for the different county sites ( Table 5). The final saturation point/volumetric soil water content (SAT) values ranged between 0.449 cm 3 /cm 3 and 0.690 cm 3 /cm 3 (Table 5). These increases in delta (DUL-LL) and SAT were done with reduction in SLPF that led to low RMSE percentages and greater than 0.5 d-statistic. In all the sites, Soil Root Growth Factor (SRGF) varied throughout the layers ( Table 6). The slightly high RMSE percentage indicated an uncertainty in soybean annual yield simulations, possibly associated with biotic stresses that the model does not simulate (disease, insects, etc.). These uncertainties are reflections of factors and events that are not considered in the model (Liu et al., 2011). The final calibrated SDUL and SLL values were within the ranges for texture and water retention data results by Fatliff et al. (1983). The use of generic DSSAT soybean genetic coefficients in this study together with the other assumptions may have contributed to the higher uncertainties.
Model R square values was ignored during the assessment of model performance because the study used time series data (30-year historic yields) which are auto-correlated values. To overcome this, correlation-based statistics (d-statistic and RMSE efficiency measures) were used to assess model performance. The yield simulations closely mimicked observed yields for the De Kalb and Dallas Counties while for Limestone and Baldwin counties the model slightly overestimated during the high yielding years and underestimated for low yielding years. The greater than 0.5 d-statistic, indicate that the CROPGRO-Soybean model yield simulations were in good agreement with observed yields.

Changes in Projected Climate Variables
To assess the changes in projected climate variables, the baseline period  representing current climatic conditions was used as a reference for comparison  for all four selected counties ( Table 7). The future annual average maximum and minimum temperatures and precipitation for the eight (8) IPSL-CM5A-MR and MIROC 5 scenarios were analyzed: 2045RCP 4.5, 2045RCP 8.5, 2075RCP 4.5 and 2075. All four sites were projected to experience increases in both maximum and minimum average annual temperatures between 1˚C to 4˚C depending on site, scenario, and year. However, the trend of change in projected rainfall patterns is not uniform for all study sites. Rainfall for Limestone and De Kalb Counties is projected to increase in the future scenarios in different amounts. Dallas County is projected to experience the highest rainfall increase in 2045 RCP 4.5 of 196 mm (MIROC 5) and 117 mm (IPSL-CM5A-MR). Rainfall for Baldwin County is projected to decrease under the extreme scenarios. Climate projections by MIROC 5 appear to be less extreme when compared to projections from IPSL-CM5A-MR.

Climate Change Impact Assessment
Simulation results reflect the yield responses of the three maturity groups (MG5, MG6 and MG7) to the different projected climate scenarios and different county soils. The MG5 yields were simulated for both De Kalb and Limestone Counties while MG6 and MG7 maturity groups were simulated for Dallas and Baldwin Counties, respectively. Rain-fed production without CO 2 for the years 2045 and 2075 was first analyzed followed by an analysis of rain fed production with CO 2 fertilization for the year 2075 when carbon dioxide concentration was projected to exceed ambient concentrations by 100 ppm. Percent changes in projected av-erage future soybean yields in comparison to current baseline average yields are shown in Figure 4 and Figure 5. The simulations for the year 2045 show soybean yield decreases in all counties for both climate scenarios (Figure 4). Projected yields decreased between 6% to 64% and between 4% to 41% for 2045 RCP 4.5 and RCP 8.5 respectively. The simulations for the year 2075 show soybean yield decrease over all the counties in both climate scenarios ( Figure 5). Yields is projected to decrease between 5% to 43% in RCP 4.5 and between 16% to 64% in RCP 8.5.
When CO 2 fertilization was introduced into the IPSL-CM5A-MR scenario soybean yield simulations of 2075, the decreases were lower compared to the simulations without CO 2 (Figure 6). These decreases ranged approximately 12% to 19% lower in the medium emission scenarios (RCP 4.5) and between 8% to 18% lower in the high emission scenarios (RCP 8.5) when compared to simulations without CO 2 fertilization. The CO 2 fertilization led to yield increases of 10% in Limestone County in the medium emission scenario. Overall, the CO 2   fertilization effect resulted in projected average yield decreases of between 9% to 41% under RCP 4.5 and RCP 8.5 respectively.
When CO 2 fertilization was introduced into the MIROC 5 scenario soybean yield simulations of 2075, the decreases were also lower compared to the simulations without CO 2 (Figure 7). These decreases were approximately 15% -23% lower in the medium emission scenarios (RCP 4.5) and 22% -32% lower in the high emission scenarios (RCP 8.5) when compared to the "without CO 2 fertilization" simulations.
The CO 2 fertilization led to yield increases of 14% and 12% in Limestone County in RCP 4.5 and RCP 8.5 respectively. Overall, the CO 2 fertilization effect resulted in projected average yield decreases of 7% and 9% under RCP 4.5 (medium emissions) and RCP 8.5 (high emissions) respectively for the time period 2075 (2060-2090). In Limestone County, carbon dioxide fertilization effect was high enough to compensate negative effects of climate change resulting to projected yield increases of 14% and 12% in the medium and high emission scenarios, respectively. These results are in agreement with other studies which show that increasing atmospheric CO 2 could have a compensation effect on yield decreases (Southworth et al., 2002;Res, Brassard, & Singh, 2007).

Discussion
The research results highlight the unique response of soybeans to increasing temperatures and the model's ability to capture these responses (Tsuji et al., 1998). In soybeans, the vegetative and reproductive stages co-exist for some period of the crop life cycle with their vegetative growing periods lengthening once temperatures exceed their optimum i.e. 22˚C -28˚C (Boote, 2011) and the crop is exposed to these temperatures continuously (Setiyono et al., 2007). Observations from field experiments, showed that this temperature-induced lengthening of the life cycle leads to an improvement in sources (Kumagai & Sameshima, 2014), i.e. increase in leaf area and leaf photosynthesis which contributes to an increase in photosynthesis and assimilate partitioning which results in increased yields. However, Thuzar et al. (2010) showed that the same increase in temperature that prolongs and enhances soybean vegetative growth and development also depresses reproductive growth leading to a decrease in the number of seeds produced. Findings by Setiyono et al. (2007) showed that reproductive growth would only be depressed if increased temperatures are experienced during flowering as there was no correlation between mean temperature and post flowering duration. The response of the different developmental stages to temperature and the influence of the timing when increased temperatures are experienced all combine to explain the inconsistent response in the medium and high emission scenarios of the soybean maturity groups simulated in Limestone and Baldwin counties. Additionally, cultivar specific day-length requirements influence the manner in which different maturity groups of soybeans respond to increasing temperatures (Kumagai & Sameshima, 2014). Boote (2011) highlighted this relationship when he showed that for early maturing groups, increasing temperatures would lead to a shorter vegetative growth period, smaller leaf area, earlier flowering and pod set and eventual decreased yields. Additionally, his simulations showed that yield reductions by increasing temperatures would be smaller in late-maturing soybean groups (higher maturity groups and higher day length sensitive groups), creating an impression that late-maturing soybean cultivars would benefit from projected future temperature increases. This explains why MG7 simulated in Baldwin County on soils with poor water holding capacity still showed resilience under climate change.
The CROPGRO-Soybean model's soybean phenological simulations show how increased temperatures in future climate scenarios affected both the time of anthesis (and subsequent seed formation) and the time to maturity (Figure 8 and Figure 9). Both shortening and lengthening of the growth stages were experienced depending on whether the temperatures exceeded the optimum soybean temperatures (22˚C -28˚C). The anthesis period was lengthened by 1 to 6 days depending on site location, GCM, RCP scenario and simulation year. The shortened development stages resulted in shortened life cycle, which means that crops will have less time to make use of available resources resulting in decreased yields. Similar studies found that hastened maturity and shortened life cycles contribute to decrease in yields (Tubiello et al., 2002;Southworth et al., 2000;Wang et al., 2011). Lengthening of the development stages, especially the time to anthesis, contributes to increase in yields or less yield decreases observed. As the vegetative and reproductive stages co-exist, the overall time to maturity was not greatly changed.  The CROPGRO-Soybean model simulations revealed the different inter-and intra-specific responses of different soybean maturity groups to climate change. MG5 simulated in Limestone County with Decatur silt loam, a soil with a high water holding capacity, showed less vulnerability to climate change, while MG7 simulated in Baldwin County with Dothan sandy loam, a soil with a low water holding capacity, also showed resilience to climate change as it is a late maturing cultivar.
This study employed a stepwise hindcast method to estimate spatially variable soil parameters in calibrating the CROPGRO-Soybean model. The model was run to assess its performance with different combinations of soil parameters and a 1:1 scatter plot was used to determine and select the set of parameters that gave the highest d-statistic and lowest Root Mean Square Error (RMSE). The study found that estimating soil parameters using aggregated regional data (in this casecounty level data) resulted in calibrated model simulations that satisfactorily mimicked long-term regional yields under current climate conditions, thus ensuring reliability in regional climate change studies. This method was practical for calibrating crop models to historical yields using aggregated regional data sets.
The calibrated crop models showed that climate change will adversely affect soybean production under both medium (RCP 4.5) and high (RCP 8.5) radiative forcing conditions used in the study. Soybean yield was projected to decrease by an average of 29% and23% in 2045, and19% and43% in 2075, under RCP 4.5 andRCP 8.5, respectively. The projected decreases were lower in MIROC 5 scenarios compared to IPSL-CM5A-MR scenarios, which had extreme projections of increase in temperature and decrease in rainfall. Considering current ongoing global efforts to mitigate climate change, the study concludes that realistic climate change impacts on soybean yields in Alabama would more likely to be within projections obtained under RCP 4.5 scenarios. Phenological simulations show that yield decreases were mainly caused by shortening crop life cycles. However, certain factors counteracted the impact of shortened life cycles on crop yields. Factors such as soil water holding capacity, soybean maturity group, timing of above optimum temperatures and water availability during the growing season, influence the crop response to climate change. The overall negative yield projection is a clear call for the development of adaptation strategies and policies to sustain current crop yields. This may involve the adaptation of new maturity groups, and development of improved cultivars, which are heat tolerant.