Evaluating the Agricultural Production Systems sIMulator (APSIM) wheat module for California

Field data

Both public and internal data including grain yield, grain protein and phenology for common wheat (Triticum aestivum L.) were obtained from the state-wide genotype trial results of the UC Small Grains Program, USA (George et al. 2017; Nelsen et al. 2018, 2019). These comprised 12 field locations in the main cereal production regions of California, conducted over three winter seasons from 2016–17 to 2018–19 (Fig. 1, Table 1). In total, 34 environments (location-by-management-by-year combinations) were sampled. The field locations fall between latitudes 33°N and 42°N, a north–south distance of ~800 km. Two locations in the dataset were identified where field yields were considered unlikely in view of the reported environmental and management conditions for the locations as described below. Estimates of the water-limited yield potential, using Sadras (2020), suggested the reported yield at Tulare in 2017 was too high given the reported growing-season rainfall and starting soil water. The yields at the low-nitrogen managed-stress Fresno 2018 site were also considered too high considering the reported nitrogen fertilisation and residual soil nitrogen. Owing to the uncertainty around the reliability of reported environmental and management variables at these locations, they were removed from the simulation.

Fig. 1.

Locations used as a source of field data to test the performance of the APSIM-Wheat module in California (George et al. 2017; Nelsen et al. 2018, 2019).

Environmental and management details are summarised in Tables 1, 2 and 3. Additional trial methodology and field management details are summarised by Nelsen et al. (2019) and George and Lundy (2019). A range of management practices were used across the field sites, including low-input dryland production, with low growing-season rainfall, to high-input and fully irrigated production. At the Davis and Fresno locations, managed-stress trials were also conducted adjacent to conventionally managed trials with the same wheat genotypes. The managed-stress trials consisted of either restricted irrigation or no nitrogen fertilisation.

Soil types at the field locations were predominantly loams, with relatively undifferentiated profiles to a depth of at least 1 m (California Soil Resources Lab 2020). Physical and chemical information for the soils at the field locations was obtained from soil samples taken during the establishment of field sites and the previous work of George and Kaffka (2017). Additional soil information was obtained from the California Soil Resources Lab (2020). Representative starting soil-water data were obtained from soil sampling and reported field measurements (George et al. 2017; Nelsen et al. 2018, 2019). The plant-available soil-water content was estimated by using this information and the crop and soil lower limits. In most locations, starting soil-water content was small relative to rainfall and irrigation.

Weather data were obtained from the California Irrigation Management Information System (CIMIS) weather station network, using the nearest representative station to each location (CIMIS 2020). The locations represent a range of climate types (Mediterranean to desert) (Peel et al. 2007). Growing-season rainfall was estimated from sowing to 1 July of the following year. Growing-season rainfall varied from 0 to 800 mm. Growing degree-days were estimated with daily average temperature, using base 0°C.

Twenty commonly grown wheat genotypes were tested across all three growing seasons (Table 4); however, not all genotypes survived to harvest in all environments. The dataset comprised a range of maturity types and yield potentials, with field yields ranging from 22 to 11 700 kg/ha.

At the Davis location in the 2016–17 and 2017–18 seasons, phenological observations were taken of the genotypes SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657). Commencing at early tillering, observations were taken at approximately ten-day intervals from plots of both the conventional and managed-stress trials. In the 2017–18 season, observations of heading and anthesis were also taken at approximately weekly intervals from all genotypes at the Davis field location. Cumulative thermal time (degree-days) spent in individual phenological stages was estimated for SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657) at the Davis location in the 2016–17 and 2017–18 seasons.

APSIM procedures

Wheat production was simulated using the unmodified APSIM-Wheat module (ver. 7.10 r4220) (Holzworth et al. 2014). Model parameterisation followed the methods for canola in California, previously described by George and Kaffka (2017), using data for the sites presented in Tables 1, 2 and 3. The field data, and previous work of George and Kaffka (2017), were used to parameterise the APSIM soil module to reflect the properties of the soil at the test locations. Sites were all sown on a fixed date according to the reported sowing date for the location in Table 1, with sowing parameters matching those described by George and Lundy (2019). Locations were fertilised and irrigated according to the reported field management in Tables 2 and 3. The APSIM Climate Control module was used to specify an atmospheric carbon dioxide content of 410 ppm, representing the approximate atmospheric carbon dioxide content for the 2016–17, 2017–18 and 2018–19 period (Blunden and Boyer 2021). To compare field observations with APSIM predictions, grain production and grain protein at harvest were simulated for all environments using cultivars that represent all unique combinations of cultivar parameters in the APSIM-Wheat module. The phenological development of all cultivars was simulated on a daily basis for the 2016–17 and 2017–18 seasons at the Davis location for comparisons with field observations of phenology taken at this location during these growing seasons.

Genetic parameters used to simulate the cultivars in the APSIM-Wheat model referred to in this work, and a description of their meaning, are presented in Tables 5 and 6. For additional information regarding model parameters, please refer to the documentation for APSIM-Wheat Module (ver. 7.10) (Zheng et al. 2015).

Data management and analysis

Data management and analyses were performed using the R software (v3.6.2) (R Core Team 2020). The accuracy of APSIM simulations for simulating grain yield and protein content was assessed by individually comparing the yields predicted by APSIM cultivars with the reported yield of wheat genotypes tested in the state-wide trials. The nature of the linear relationship between the simulated and field data was assessed using the root-mean-square error (RMSE), the coefficient of determination (R²), and the slope and intercept of the linear relationship (Ahmed et al. 2016; Brown et al. 2018; Wallach et al. 2019). For comparison with published literature, the RMSE was also normalised (nRMSE) using the data range method of the hydroGOF library in R (Zambrano-Bigiarini 2020).

Management and environmental variables were unique to individual locations, so their impact on model output was explored by comparing the difference between simulated and observed yields at the location level with management and environmental variables across locations. The relative importance of individual APSIM-Wheat module genetic parameters (specified in Table 5) for the performance of the simulation was assessed by using a classification and regression tree from the rpart library in R (Therneau and Atkinson 2019), and by stepwise linear regression, implemented using the olsrr library in R (Hebbali 2020). For classification and regression tree (CART) analysis, genetic parameters for individual APSIM cultivars were used as predictor variables of the RMSE between the simulated and actual yield and protein data across all environments. A grid search approach was used to determine optimal model parameter settings for the CART analysis that minimised the estimates of cross-validated prediction error. Stepwise linear regression was performed assuming genetic parameters were numeric variables. Given differences in units, the data were scaled and centred. The best model from the stepwise linear regression was selected based on the Akaike Information Criterion (AIC). The stability of regression estimates was assessed using variance inflation factors, which measure the inflation in the variances of the parameter estimates due to collinearities that exist among the predictors (Hebbali 2020). The regression diagnostics of initial stepwise linear regression models found that node senescence on the main stem exhibited multiple collinearities with other predictors, as well as residual-leverage issues, and it was therefore removed from the stepwise regression analysis. Initial analyses also found that interactions between parameters could not be included in the regression model owing to collinearities and model sparsity. Regression assumptions of the linear fits were visually assessed. The process described above was also applied to the coefficient of determination; it did not lead to different conclusions from the RMSE and is therefore not reported.

Model calibration

Field data needed for comprehensive model calibration and validation were not available. However, using findings from the analyses described above, preliminary cultivar parameterisation was performed. The parameter file for wheat was modified to match the thermal-time values observed for SY Cal Rojo (1478) and SY Blanca Grande 515 (1657) reported in Table 6. This was done for the base cultivar as well as the New Zealand base cultivar, which generated some of the more reliable simulations and these are described in the Results.

Performance of the APSIM-Wheat model for simulation of yield and protein in California

The accuracy and precision of the existing APSIM-Wheat module for predicting grain yield and protein varied considerably between unique combinations of APSIM-Wheat cultivars and field genotypes (Fig. 2). For grain yield, across all pair-wise comparisons, the relationships observed for RMSE ranged from 1040 to 4900 kg/ha and nRMSE from 16% to 68%. The RMSE of protein varied from 1.4% to 5.1% and the nRMSE varied from 17% to 74%.

Fig. 2.

Distribution of root-mean-square-errors (RMSE) between all combinations of APSIM-Wheat cultivars and field genotypes tested in the analysis for (a) grain yield and (b) grain protein. Dashed line shows the mean value.

In terms of grain yield, the best linear relationships between the simulated and field data for SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657) are presented in Fig. 3a, b. In both cases, the nRMSE was >20%, with multiple outliers, and the linear fit over-predicted lower yielding environments and under-predicted high-yielding environments. Across all pairwise combinations of APSIM cultivars and field genotypes, WB Joaquin Oro (UC1728) had one of the lowest values of RMSE for yield (Fig. 3c). There were few apparent outliers, and the fit was close to 1:1, despite over-prediction in lower yielding environments. The genotype LCS Atomo (UC1723) had the highest yields of all genotypes tested in the multi-environment trials and was therefore specifically examined. The best model fit for this genotype was relatively poor (Fig. 3d). The relationships between simulated and field data for grain protein content were considered poor by all metrics, the best fit being r² = 0.35, and are not presented.

Fig. 3.

Relationships between grain yields from field and simulated genotypes: (a) SY Cal Rojo (UC1478) and APSIM cultivar Gatton Hartog (RMSE = 1750, nRMSE = 22%, mean absolute error = 1300, r² = 0.30, intercept = 2720, slope = 0.48); (b) SY Blanca Grande 515 (UC1657) and APSIM cultivar Gregory (RMSE = 1730, nRMSE = 21%, mean absolute error = 1346, r² = 0.32, intercept = 2512, slope = 0.48); (c) WB Joaquin Oro (UC1728) and APSIM cultivar V2_P2 (RMSE = 1129, nRMSE = 17%, mean absolute error = 1000, r² = 0.59, intercept = 1257, slope = 0.77); (d) LCS Atomo (UC1723) and Gatton Hartog (RMSE = 2084, nRMSE = 24%, mean absolute error = 1694, r² = 0.27, intercept = 2599, slope = 0.44). Solid line shows the unconstrained linear relationship and dashed line the 1:1 relationship. Shading indicates the 95% confidence interval.

Performance of the APSIM-Wheat module for simulation of phenology

Figs 4 and 5 show the relationships between the observed field phenology of SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657), respectively, at the Davis location, and APSIM model predictions of phenology using the cultivars that gave the smallest RMSE of prediction for yield (as shown in Fig. 3). In all environments, the predicted time to all phenological stages between tillering and anthesis was slightly shorter than observed in the field, and was longer for late anthesis, early grain development, and grain maturation. The phenological data for WB Joaquin Oro (UC1728) and LCS Atomo (UC1723) showed a similar pattern but were available only for the Davis location in a single season, and at a lower temporal resolution, and are therefore not presented.

Fig. 4.

Relationship between observed phenology of field genotype SY Cal Rojo (UC1478) and predicted phenology of simulated APSIM cultivar Gatton Hartog (filled circles). Field observations span a range of phenology stages and are therefore represented as a value range (vertical lines). Observations were taken at the Davis location in the 2016–17 and 2017–18 seasons. ln, low-nitrogen managed-stress trial; lw, low-water managed-stress trials.

Fig. 5.

Relationship between observed phenology of field genotype SY Blanca Grande 515 (UC1657) and predicted phenology of APSIM cultivar Gregory (black circles). Field observations span a range of phenology stages and are therefore represented as a value range (vertical lines). Observations were taken at the Davis location in the 2016–17 and 2017–18 seasons. ln, low-nitrogen managed-stress trial; lw, low-water managed-stress trials.

Importance of management, environmental, and APSIM-Wheat parameters to simulation reliability

Water and nitrogen status of the locations was not strongly associated with model reliability; there was a trend whereby the yield at sites with lower water and nitrogen status was under-predicted, although the relationship was weak with large residual variation (Figs 6 and 7). CART analyses found thermal time to the start of grainfill and number of grains per stem weight were strong determinants of model reliability, as were vernalisation and photoperiod sensitivity (Fig. 8). Stepwise linear regression found all variables to have a significant impact on RMSE (Table 7), and similar to the CART analysis, thermal time to the start of grainfill and number of grains per stem weight had some of the largest impacts on model reliability. By contrast, vernalisation and photoperiod were not found to be strong predictors. The analysis was also performed using the coefficient of determination as the indicator of model reliability, but the results were the same as for RMSE between simulated and actual yield and protein data, and the results of the analysis are therefore not presented.

Fig. 6.

Relationship between total water (starting plant-available soil water, seasonal precipitation and irrigation) and the difference between predicted and observed yield across all locations and years. Dotted line shows no correlation between predictor and response (y = 0).

Fig. 7.

Relationship between total nitrogen (starting soil nitrogen and fertilisation) and the difference between predicted and observed yield across all locations and years. Dotted line shows no correlation between predictor and response (y = 0).

Fig. 8.

Classification and regression tree showing the importance of predictor variables to the RMSE of the linear relationship between all combinations of APSIM wheat cultivars and field genotypes for grain yield. See Tables 5 and 6 for a description of the values. Boxes along the bottom depict RMSE and the proportion of the dataset that falls in each division.

Observations of thermal time for phenological stages and model calibration

The thermal time spent in phenological stages was similar for both SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657), under conventional and managed-stress conditions in individual seasons. Comparison of thermal time spent in phenological stages with the equivalent parameter in the APSIM-Wheat module (Table 6) shows that the thermal times needed from sowing to the end of juvenile stage, from floral initiation to flowering, and to reach the anthesis stage were greater than for the APSIM cultivars, whereas thermal time needed from the beginning to end of grainfill was smaller.

Cultivar calibration

The APSIM base cultivar and New Zealand base cultivar were both reparameterised with thermal-time observations from SY Cal Rojo (UC1478) and SY Blanca Grande 515 (UC1657). However, like the unparameterised cultivars reported in Fig. 3, the reparameterised cultivars simulated WB Joaquin Oro (UC1728) most reliably (Table 8). The summary statistics show that the relationship between predicted and observed yield for the parameterised cultivars was worse than for the unparameterised model in Fig. 3.

The performance of APSIM for simulation of wheat yield in California

Computer-based crop simulation models support agricultural research and production, but model calibration and testing requires field data that are often not available. The objective of this work was to use pre-existing data from state-wide field trials to evaluate the performance of the APSIM model (v7.10) for simulating wheat production in California and to identify future research needs for model validation and calibration. Given the large yield range and the diversity of environments, management methods and latitudes represented in the field data, we found that the APSIM model was capable of predicting relative grain yield responses for at least one field genotype. The accuracy of the relationship between predicted and reported field yields was comparable to other tests of the APSIM-Wheat module in environments where it has not been previously calibrated (Asseng et al. 2000; Farré et al. 2002; Zhao et al. 2014, 2020; Kouadio et al. 2015; Hussain et al. 2018).

Despite this, we do not consider the current APSIM-Wheat module to be accurate or precise enough to simulate wheat production scenarios in this environment, where greater specificity is required. The model could not reliably simulate high-yielding genotypes and was unable to reliably simulate grain protein content. The genotype that was simulated most reliably, WB Joaquin Oro (UC1728), was one of the lowest yielding of the advanced genotypes included in the California field trials, although it achieved some of the highest protein values (Nelsen et al. 2018, 2019). The usefulness of the current APSIM model for the Californian wheat industry is therefore limited, and it will require further calibration and validation.

Reasons for the performance of APSIM-wheat in California

The poor performance of APSIM for simulating wheat production in California is surprising given the environmental similarities between the region and southern Australia, where the model has been calibrated and extensively validated. The uncalibrated APSIM model has also previously been found to simulate canola production reliably in the same region (George and Kaffka 2017). Inaccuracies in crop simulation modelling result from both parameter sensitivity and parameterisation error (He et al. 2015).

Correct calibration of APSIM relies on genetic parameters that accurately represent the crop genotype being simulated (Ceglar et al. 2011; Zhao et al. 2014; Casadebaig et al. 2016; Harris et al. 2016; Meier et al. 2020). As would be expected, the performance of APSIM for simulating wheat in California was shown to vary considerably depending on the choice of APSIM cultivar and the field genotype with which it was compared. Analyses of the importance of variables for model reliability, observation of field phenology, and initial thermal-time estimates all suggest that the parameterisation of APSIM cultivars does not currently represent the Californian genotypes well. APSIM was able to approximate the phenology of field genotypes SY Cal Rojo (UC1478), SY Blanca Grande 515 (UC1657) and LCS Atomo (UC 1723) despite the yields of these genotypes not being simulated reliably. By contrast, the phenology of WB Joaquin Oro (UC1728) was not simulated as reliably, despite the yield of this genotype being simulated more reliably than that of other genotypes.

Several phenology and grain parameters are influencing the outcome of the simulation. Parameterisation of the model using field observations of thermal time did not lead to improvements in yield simulations. Therefore, fieldwork is necessary to develop a better understanding of these phenology and grain parameters in Californian varieties. Most wheat genotypes in California are considered to have a low sensitivity to vernalisation and are photoperiod-insensitive (Dubcovsky et al. 2006). However, Ottman et al. (2013) demonstrated that consideration of photoperiod alongside thermal time improves predictions of flowering date in autumn/fall-planted spring wheat cultivars in the south-west of USA. Our findings also suggest that photoperiod and vernalisation sensitivity are important predictors and should also be investigated in more detail.

All efforts were made to parameterise the model to accurately reflect environmental and management conditions at the locations using reported data, but there are reasons to believe the management and environmental data for some locations are not accurate. The RMSE of prediction of yield of WB Joaquin Oro (UC 1728) was also comparable to the reported single-site standard deviation from the field data for the genotype (George et al. 2017; Nelsen et al. 2018, 2019), illustrating the sensitivity of yield to environmental factors within field sites not being captured by the simulation, along with the potential for imprecisions in reported management and environmental data at field sites.

Datasets from state-wide multi-environment trials, which are often long-running, geographically extensive, agro-environmentally diverse, and publicly available, are resources for understanding the impacts of management and environmental data on crop production (e.g. George and Lundy 2019; Barrett-Lennard et al. 2024). Within resourcing limits, the California Small Grains Program collects comprehensive data regarding crops and test locations (George et al. 2017; Nelsen et al. 2018, 2019), but our study shows such datasets may still have limitations due to insufficient resolution or completeness. Well-conducted multi-environment trials are essential for effective plant breeding and agronomy (Yan 2014) and will continue to be routine in many regions of the world. Existing multi-environment trials should therefore be leveraged to generate data that can be applied to analytics and modelling (Keating 2020; Pasquel et al. 2022). This will increase the long-term value of these programs.

Conclusions and future work: improvement of the APSIM-Wheat module for California

The current APSIM-Wheat module calibration is not able to simulate wheat production reliably under Californian conditions. Models such as APSIM will continue to be important tools for agricultural research and management (Chenu et al. 2017; Keating 2020), and increasing the accuracy and precision of the APSIM-Wheat module for Californian conditions will provide a tool to support and complement existing agricultural research. Our study finds that data needed for model calibration are not currently available. Future work should therefore prioritise higher resolution observations of phenology, sensitivity to photoperiod and vernalisation, and grain characteristics, as well as canopy development and senescence (Zheng et al. 2015; Brown et al. 2018). Seasonal biomass accumulation and soil-water dynamics would also be valuable (Asseng et al. 1998) and assist in improving simulations.