Probabilistic Assessment of Cereal Rye Cover Crop Impacts on Regional Crop Yield and Soil Carbon

Abstract

Field research for exploring the impact of winter cover crops (WCCs) integration into cropping systems is resource intensive, time-consuming and offers limited application beyond the study area. To bridge this gap, we used the APSIM model, to simulate corn (Zea mays L.)-rye (Secale cereale L.)-corn-rye and corn-rye-soybean (Glycine max L.)-rye rotations in comparison with corn-corn and corn-soybean rotations across the state of Illinois at a spatial resolution of 5 km × 5 km from 2000 to 2020 to study the impact of WCCs on soil organic carbon (SOC) dynamics and crop production. By propagating the uncertainty in model simulations associated with initial conditions, weather, soil, and management practices, we estimated the probability and the expected value of change in crop yield and SOC following WCC integration. Our results suggest that integrating cereal rye into the crop rotations imparted greater yield stability for corn across the state. It was found that the areas with low probability of increase in SOC (p < 0.75) responded equally well for soil carbon sequestration through long term adoption of WCCs. This study presents the most complete uncertainty accounting of WCC benefits across a broad region and provides greater insights into the spatiotemporal variability of WCCs benefits for increasing WCC adoption rate.

1. Introduction

Corn (Zea mays L.) and soybean (Glycine max L.) production in the U.S. Midwest is an undeniably vital component of U.S. and global food systems, accounting for more than 80% of national and 25–30% of global output for these commodities [1]. Agricultural lands are typically fallowed in the winter months in the U.S Midwest following intense summer growing seasons. This setup poses significant agronomic and environmental challenges such as water quality issues across the U.S Midwest [2,3] and the formation of the hypoxic zone in the Gulf of Mexico. Adopting conservation agricultural practices through integrated nutrient management and complex crop rotations has been reported to overcome some of the challenges with respect to nutrient recycling and enhance carbon sequestration [4,5]. Another common tool that can play a key role in soil conservation and future carbon markets is the integration of winter cover crops into current rotations [6,7]. Some of the widely reported benefits of WCCs include soil organic matter (SOM) enhancement, improved soil infiltration rates, and most importantly, reduction in nutrient losses and soil erosion [8]. The associated benefits greatly depend on the choice of WCC type, and cereal rye (Secale cereale L.) is one of the most widely used species in U.S. production systems [9]. Additional advantages of cereal rye include its extreme cold-hardiness and ability to produce considerable biomass, thus, making it a popular choice for protecting soil in intensive cropping systems in the U.S. Midwest [10].

However, WCCs have not been widely adopted, largely due to considerable uncertainties regarding their economic, agronomic, and climate benefits [11,12]. These uncertainties stem from the highly uncertain role of WCCs in altering the water and carbon dynamics of cropping systems, which then impacts the subsequent cash crop and the farmer’s main source of income. For instance, a meta-analysis did not find any statistically significant difference within the winter rye data relating WCCs with yield increases [13], while more recent studies have even shown a six percent reduction in corn yield following a winter rye cover crop [14]. Advocacy programs for increasing WCCs adoption have been developed partially due to the preliminary data suggesting an increase in soil carbon sequestration and reduction in nitrate leaching as the two most important climate and environmental services of WCCs [12,15]. However, each of these benefits is strongly dependent on the environment, as WCCs have also been shown to excessively take up water from soil in low precipitation areas, potentially exacerbating the drought condition [16]. Further, little is known about the capacity of WCCs to sequester carbon across different environments. Consequently, there are highly variable reports on the relationship between WCCs and carbon sequestration in the literature, adding an extra layer of agronomic uncertainty. The majority of reported experiments [17,18] assessing the effect of WCCs on soil organic carbon (SOC) were often conducted only for a few years, and given that annual rates of increase in SOC are greatest in the early years, it is likely that the current suggested rates of SOC increase are substantially overestimated [19]. In addition, climate and environmental services of WCCs have been shown to have strong spatial and temporal uncertainty due to a multitude of environmental conditions and associated interactions between climate and the processes that influence WCC growth and development [20].

While small-plot research experiments are essential in a variety of cropping systems across the U.S Midwest, large scale and long term exploration of climate and agronomic benefits of WCCs is especially needed for different climates and environments. The lack of regional and long term observational data on the performance of WCCs has hindered the fast-paced research required to provide necessary information to farmers for large-scale adoption [21]. Process-based crop models such as the APSIM model (Agricultural Production Systems Simulator) [22] can fill this gap by providing a system level representation of different soil and crop processes with explicit representation of crop genetics. APSIM is a modeling program with multiple biophysical modules based around, soil, weather, and management. These pioneering models simulate complex cropping systems in an internally consistent manner by conserving mass and energy [23,24] and have been extensively used in the literature to explore the complex interactions between G × E × M across the world [25]. In the past, process-based crop models have been used to assess the impact of WCCs on crop production, soil water dynamics, and greenhouse gas emissions [26,27]. However, most studies were limited to a single or few sites [1,28] and only one aspect of WCC impacts on soil and crop dynamics was explored. Consequently, most of these studies fall short of capturing the spatial variation in WCC impacts across broad regions and miss the opportunity of identifying the drivers of this variation. Furthermore, most of the previous WCC modeling studies lack proper accounting of the uncertainties associated with model inputs and parameters such as weather, soil properties, and management practices. Therefore, such studies, although beneficial for their local sites, are not reliable for understanding the system at broader scales due to lack of generalizability. Considering the ample potential that WCCs hold, we provide the first and most complete uncertainty accounting of WCCs’ potential impacts on improving crop yield and SOC sequestration across the state of Illinois. The main objectives of the current study are to: (a) constrain uncertainty in APSIM model parameters by calibrating and validating the APSIM-Wheat model for cereal rye; (b) quantify the potential for WCC biomass accumulation across the state of Illinois; and (c) assess the impact of WCC integration on SOC and cash crop performance.

2. Materials and Methods

2.1. Study Area

The study area for this crop modeling experiment was the entire state of Illinois, a key contributor to the U.S. Corn Belt (Figure 1). The average area under maize production for the study period of 2005 to 2020 was around 4.8 ± 0.3 million ha, with an average yield of 10929 ± 1623 kg ha⁻¹. Similarly, the average area under soybean production in Illinois was estimated to be around 3.9 ± 0.3 million ha, with an average yield of 3486 ± 426 kg ha⁻¹ [29]. Across the study period, 2005 was the driest year with respect to cumulative seasonal precipitation (464 mm) (April–October), and 2019 was the wettest year (935 mm). The year 2010 recorded the hottest growing season (April–October) based on average daily temperatures, whereas 2009 was found to be the coldest.

2.2. Simulation Setup

2.2.1. Modeling Platform

The current study was based on APSIM v7.9 which is one of the most popular open-source cropping system models. The model can simulate various crops across diverse geographical regions and their interactions with soil and atmospheric conditions at a point scale at daily timesteps. As part of the input data, information on crop parameters, management practices, soil properties, and weather is required for simulating crops in APSIM. The maize model in APSIM is based on the CERES-Maize model [30], whereas soybean is simulated using the generic PLANT model native to the APSIM platform (for more information on APSIM, see https://www.apsim.info/; accessed on 15 April 2021). We used APSIM-Wheat, which also uses the PLANT module, as a proxy for simulating cereal rye and will be referred to as the ‘rye model’ hereafter.

We used the parallel system for integrating impact models and sectors (pSIMS) [31] modeling platform to run APSIM simulations on a regional scale which established a uniform protocol for the aggregation and harmonization of input and output data. In addition, pSIMS enabled the parallel implementation of APSIM simulations on a large spatial grid and facilitated post-processing of the simulated data.

2.2.2. Uncertainty Propagation

To account for all sources of uncertainty (i.e., model structure, parameters, and inputs) that contribute to the overall uncertainty in model predictions, we defined five broad uncertainty classes that included initial conditions, soil, weather, management, and crop parameters (

Table 1. Uncertainty factors considered for simulating each scenario in long term simulations. U stands for uniform distribution, and N stands for normal distribution.

Name	Options	Definition
Initial Conditions	Residue type (RT)	RT ~ sample (corn, soybean)
	Residue weight (RW; kg ha−1)	RW ~ U (100, 2500)
	Water fraction (WF)	WF ~ U (0.05, 0.95)
Soil	GSDE/SoilGrid
Weather	10 ensembles from ERA5
Management	Planting date (pdate)	pdate + N (µ = 0, σ = sd(pdate))
	Harvesting date (hdate)	hdate + N (µ = 0, σ = sd(hdate))
	Rye seeding rate (plpop; seeds m−2)	plpop ~ U (200, 500)
Parameters	Corn: Ensemble of 6 cultivar parameters
	Soybean: Ensemble of predefined cultivars depending upon maturity group, based on latitude (30 total genotypes varying from MG 2 to MG 4)
	Rye: Ensemble of 7 optimized genotypes

) [32]. We included uncertainty propagation in all steps performed in this study, including the sensitivity analysis, emulator generation, and final regional simulations. In our study, pSIMS was set up to randomly choose a soil profile from the available input data products [31]. The global soil dataset for earth system modeling (GSDE) [33] and the SoilGrid dataset [34] at a spatial resolution of 250 m were used to ensemble soil properties. Weather uncertainty was incorporated through ten weather ensembles offered by the ERA5 reanalysis data product, a global gridded data product developed by the European Centre for Medium-Range Weather Forecasts [35]. The ERA5 weather data were available hourly at a spatial resolution of 30 km. Uncertainty around initial conditions was propagated for residue type (corn or soybean), residue weight (kg ha⁻¹), and water fraction (mm/mm) at the beginning of the simulation period. For management practices, we accounted for uncertainty in planting date and harvesting date, as well as crop parameters to account for genetic background variability.

2.2.3. Multi-Site and Multi-Criteria Model Calibration and Validation

Field data on cereal rye phenology and biomass were collected at 7 different sites across the U.S. Midwest from 2018 to 2019 (Figure 2). Cereal rye phenology was measured using Zadoks’ scale [36]. Field data, along with management inputs, were used to calibrate and validate the rye model.

Model calibration was performed by, first, identifying crop parameters that controlled crop phenology and biomass accumulation through a literature review [1,37,38,39] and a global sensitivity analysis (GSA). Rye phenology in APSIM is regulated by thermal time and is sensitive to photoperiodism and vernalization. Most crop models use radiation use efficiency (RUE) for predicting biomass assimilation in cereal crops [39,40]. Since cereal rye and other winter hardy crops undergo vernalization during the winter, it may be more effective to apply different values of RUE at different stages of the cereal rye lifecycle. To test this, we divided the cereal rye lifecycle at the crop stage code 4.0 (end of juvenile), which corresponds to the period when the temperature starts to rise during spring, and assigned a separate RUE parameter to each period to allow for temporal variation in this parameter. Moreover, the RUE parameter value for the APSIM-Wheat model [39] is much lower than what is likely required for cereal rye [40]. Therefore, by introducing temporal variations in RUE and performing a GSA analysis, we can see which RUE time period is more important for cereal rye calibration (Figure S1). In the sensitivity analysis, we selected a set of 9 parameters that directly controlled phenology, biomass, or both, and evaluated the sensitivity of those state variables to changes in values of the selected parameters (

Table 2. Parameters selected for sensitivity analysis and calibration of cereal rye phenology and biomass in APSIM.

Parameter	Description	Default Value	Priors for Calibration
pesw_germ *	Plant extractable soil water in seedling layer inadequate for germination (mm/mm)	0.00	~N (µ = 0.15, σ = 1)
tt_end_of_juvenile *	The potential period from end of juvenile stage to terminal spikelet stage (°Cd)	400	~N (µ = 400, σ = 25)
tt_floral_initiation	The potential period from floral initiation flowering stage (°Cd)	555
vern_sens *	Vernalization sensitivity	1.5	~N (µ = 5, σ = 2)
photop-sens *	Photoperiod sensitivity	3.0	~N (µ = 5, σ = 2)
y_rue1 (fall)	Radiation use efficiency for fall (g MJ−1)	1.24
y_rue2 (spring) *	Radiation use efficiency for spring (g MJ−1)	1.24	~N (µ = 2.98, σ = 1)
x_ave_temp1	Lower bound of the mean daily temperature where photosynthesis is not hindered (°C)	10
x_ave_temp2	Upper bound of the mean daily temperature where photosynthesis is not hindered (°C)	25

Note: * represents the parameter selected for calibration after sensitivity analysis.

).

An ANOVA-based GSA was conducted by taking 900 random samples from the parameter space, running the APSIM model over the 900 points across all sites, and decomposing the variability explained by each parameter to determine the most influential parameters [41]. A sensitivity index was then calculated for each parameter based on the sum of squared errors (SSQ) as follows:

$$Main effect sensitivity indices: S_1=\frac{SS{Q}_1}{SS{Q}_T}; S_2=\frac{SS{Q}_2}{SS{Q}_T}$$

$$Interaction sensitivity indices:S_{12}=\frac{SS{Q}_{12}}{SS{Q}_T}$$

$$Total sensitivity indices: T{S}_1=\frac{SS{Q}_1+SS{Q}_{12}}{SS{Q}_T}$$

where $SS{Q}_k$ is the sum of squares associated with the factor $k$, and $SS{Q}_T$ is the total sum of squares for a given variable.

Rather than using APSIM to optimize the most influential parameters, we used statistical surrogate models (emulators) to replace APSIM. Emulators are fast statistical models that can closely replicate process-based crop models in a more constrained inference space while also providing flexibility for performing more computationally expensive optimization schemes, such as Bayesian approaches [42].

After identifying the most influential parameters, 500 knots were randomly selected from the new parameter space to generate emulators for each of the seven sites. APSIM was run at each knot, and model outputs of phenology and biomass were used to build generalized additive models (GAMs) as site-level emulators to simulate biomass and phenology as a function of the most sensitive parameters. Emulators were developed in R statistical software [43] using the ‘mcgv’ package [44].

The emulators were then used in a multicriteria global Bayesian optimization scheme to constrain the parameters using all observations across all sites. Model calibration and validation was performed using a leave-one-out cross validation scheme, where a model was trained for cereal rye biomass and phenology simultaneously at six sites and tested on the outputs of the remaining seventh site. Posterior densities of model parameters were estimated via the Markov Chain Monte Carlo (MCMC) method with 5000 iterations and 500 burn-in iterations using the R package NIMBLE [45,46]. The Bayesian models were defined as:

$$\mu ~ N \left({\mu }_0, \tau \right)$$

$$Y_k^p ~ N(f_p\left(\mu ),\sigma \right)$$

$$Y_k^b ~ N\left(f_b\left(\mu \right), \sigma \right)$$

where ${\mu }_0$ is a vector of mean model parameter with $\tau$ capturing the variation from true parameters, $Y_k^p \mathrm{and} Y_k^b$ are, respectively, the phenology and biomass observation collected at different times and across different sites that are a function of the mean model parameter value, and $f()$ is the emulator. The priors were defined based on reported values in the literature (Table 2).

To compare model performance, we compared the predicted values of cereal rye Zadok stage and biomass with the observed data by estimating index of agreement (d), mean error (ME) and normalized root mean square error (nRMSE) [47], as:

$$d=1-\frac{{{\displaystyle \sum }}^{ }{\left(y_i-x_i\right)}^2}{{{\displaystyle \sum }}^{ }{\left(\left|y_i-\overline{x}|+|x_i-\overline{y}\right|\right)}^2}$$

$$ME=\frac{{{\displaystyle \sum }}^{ }\left(y_i-x_i\right)}{n}$$

$$nRMSE=\frac{RMSE}{\overline{x}}\times 100$$

where $y_i$ is the predicted value, $x_i$ is the observed value, $n$ is the number of observations, and $\overline{x}$ is the mean of observed data. The aim of the calibration was to maximize the likelihood of observing both phenology and biomass at testing sites simultaneously, which we expected to yield a higher value of d, while reducing ME and nRMSE.

2.3. Long Term Simulations

The goal of the long term model runs was to simulate the potential impact of WCCs on crop performance and SOC when integrated into different crop rotations. The modeling platform can operate at much finer or coarser spatial scales, but there is a tradeoff between the spatial resolution and computation time. To avoid the computational constraints of finer scales, the spatial resolution of this modeling study was set at 5 km × 5 km across the whole cropland area in the state of Illinois. Simulations were performed from 2000 to 2020, where the period from 2000 to 2004 served as a spin-up period to allow for the model to reach equilibrium. All results presented in later sections will focus on simulations from 2005 to 2020. Since the goal was to explore the potential benefits of a cereal rye cover crop across the region, we did not simulate the exact historic crop rotations at each pixel. Instead, we applied 4 different rotations across the study area, simulating continuous corn (CC) and corn–soybean (CS) rotations with and without WCCs. The crop rotations that included WCCs were corn–cereal rye-corn-cereal rye (CRCR), and corn-cereal rye-soybean-cereal rye (CRSR).

There exists significant uncertainty in the WCC management decisions, as cereal rye is expected to have a negative impact on corn yield and no impact on soybean yield [48]. Hence, in our simulation study, cereal rye preceding corn was terminated 14 days prior to cash crop planting, and cereal rye preceding soybean was terminated one day prior to cash crop planting. Hence, we simulated 4 different scenarios of crop rotation for this study. Each scenario consisted of 50 ensemble members to account for the uncertainty around various factors as described in Table 1. Cover crop planting was carried out one day after cash crop harvesting. Therefore, any uncertainty in the harvesting date of cash crops automatically adds to the uncertainty around WCC planting date. For crop management operations, planting dates for corn and soybean were obtained from the Agricultural Model Intercomparison and Improvement Project (AGMIP) data product [49]. The amount of nitrogen (N) fertilizer added to corn was also obtained from the AGMIP gridded data product, whereas no fertilizer was applied to soybean. The crop parameters controlling corn growth and development were taken from [50]. To account for the soybean maturity group (MG) gradient that varies from north to south, soybean cultivars were ensembled for each pixel depending on the latitude [51]. Pixels north of 41.43° N were planted with ensembles of MG 2 cultivars defined in the APSIM model, and all pixels south of 39.27° N were planted with ensembles of MG 4 cultivars. Any pixels between the two mentioned latitudes were planted with MG 3 ensembles.

2.4. Statistical Analysis

To compare the impact of WCCs on SOC and crop performance simulated by the model, we estimated the difference between ensemble means of CC with CRCR, and CS with CRSR across all pixels and years. To incorporate the uncertainty estimated around mean predictions, we calculated the probability and expected value of change in SOC and yield across different scenarios. With respect to SOC, assuming the mean difference followed a normal distribution, we estimated the probability of observing a mean difference to be greater than 0 (pSOC). Based on the estimated probability, we calculated the expected change in SOC for each pixel across all the years as:

$$ESOC=\frac{{∆}_{SOC}\times pSOC}{SO{C}_{ctrl}}\times 100$$

where $ESOC$ is the expected percentage change in SOC due to the addition of WCCs to the rotation, ${∆}_{SOC}$ is the difference between means of control and treatment ensembles, and $SO{C}_{ctrl}$ is the mean SOC for control (CC or CS rotation). The estimated pSOC and ESOC were then divided into different categories, high and low probability, and expected value. Through classification, we generated a grid of 4 classes, namely: (i) low pSOC, low ESOC (Class 1); (ii) low pSOC, high ESOC (Class 2); (iii) high pSOC, low ESOC (Class 3); and (iv) high pSOC, high ESOC (Class 4), that aided in geographical mapping of the region. The same process of estimating probabilities and expected change was repeated for corn and soybean yields to estimate probability for the increase in crop yield to be greater than 0 when WCC is integrated into the rotation (pY), and the expected change in crop yield (EY). Next, a yield stability index (YSI) was estimated for crop performance by taking the ratio of coefficient of variation (CV) for control to the CV for treatment, as

$$CV=\frac{\sigma }{\mu }$$

$$YSI=\frac{C{V}_{ctrl}}{C{V}_{trt}}$$

where $\sigma$ is the standard deviation, and $\mu$ is the mean crop yield. The YSI greater than 1 can be interpreted as an expression of greater variation in the control treatment, and vice versa. The areas with YSI greater than 1 depict regions where WCCs have imparted greater stability to the crop yield.

3. Result

3.1. Sensitivity Analysis, Model Calibration and Validation

The results of our GSA identified ‘pesw_germ’ (i.e., the moisture requirement for germination) was the most sensitive parameter, explaining around 19% of the variability in rye biomass and 20% of the variability in rye phenology (Figure S2). The next most important parameter, which explained about 8% of the variability in rye biomass was ‘y_rue2’, the RUE parameter value for spring. For variability on cereal rye phenology, other important parameters included ‘tt_end_of_juvenile’ (thermal time from end of juvenile stage to terminal spikelet stage), ‘photop_sens’ (photoperiod sensitivity), and ‘vern_sens’ (vernalization sensitivity). Our GSA analysis showed a relatively large leftover variation (Residuals in Figure S2) that was not explained by any of the examined parameters, representing significant interaction among all other factors and model parameters explaining rye phenology and biomass. This points to the need for additional observational data for constraining leftover variation.

Based on the results of the GSA, the rye model was calibrated by controlling both phenology and biomass simultaneously through optimizing the five most sensitive parameters, governing rye phenology and biomass accumulation. Crop parameters ‘pesw_germ’, ‘photop_vern’, and ‘vern_sens’ have previously been used for calibration of cereal rye biomass and phenology [38]. The resulting posterior distributions for the calibrated parameters are presented in Figure S3. Based on the comparison between the observed and predicted cereal rye data, the calibrated rye model performed well in capturing the phenology and biomass accumulation across different sites. A reasonable agreement between observed and predicted cereal rye phenology was observed with d = 0.66, ME = −5.03, and nRMSE = 17% (Figure 3). Similarly, predicted and observed biomass values showed good agreement with d = 0.84, ME = 338 kg ha⁻¹ and nRMSE = 39%. The results of validation highlight a slight underprediction of growth stages, and an overprediction of biomass. Based on the estimated validation statistics, it can be stated that the rye model performed reasonably well in capturing the growth stages and biomass accumulation across the region and can be used for scenario analysis.

3.2. Rye Biomass

For the CRCR rotation, cereal rye produced, on average, 11,793 kg ha⁻¹ of dry biomass over the study period. The lowest biomass was observed in the year 2005 (6987 kg ha⁻¹), and the highest biomass was observed in 2012 (19,206 kg ha⁻¹) (Figure S4). Greater variation in cereal rye biomass was found for the CRSR rotation. Cereal rye preceding corn (CrC), on average, produced 4806 kg ha⁻¹ in dry biomass with the highest accumulation in 2012 (6830 kg ha⁻¹) and the lowest in 2010 (2949 kg ha⁻¹). Cereal rye preceding soybean (CrS) produced greater biomass than CrC. The average biomass produced for CrS was found to be 17,878 kg ha⁻¹, with the highest biomass simulated in 2017 (24,115 kg ha⁻¹) and the lowest in 2005 (12,584 kg ha⁻¹).

We attempted to partition the rye biomass variability by considering variation in weather (seasonal precipitation, average daily temperature) and soil variables (soil texture, cation exchange capacity, pH). The most important factor which explained ~75% of the observed variation in simulated cereal rye biomass was crop rotation, followed by latitude (17%). In comparison with the rye biomass produced in CRCR rotation, CRSR produced significantly lower rye biomass when rye preceded corn as indicated by a negative slope. The other major factors contributing to the variation in cereal rye biomass were cereal rye planting and termination dates, which were directly influenced by the cash crop planting and harvesting dates. Since soybean was planted later than corn and WCCs were terminated only 1 day prior to soybean planting, CrS had a longer growing period than CrC and, thus, more time to accumulate biomass. The CrC biomass in CRSR rotation was found to be lower than rye biomass in the CRCR rotation for the corresponding years. The primary factor that contributed to lower biomass accumulation in CrC was that soybean harvesting dates showed greater variation across the state, which increased uncertainty around planting dates for CrC. Hence the growing window for CrC was greatly reduced, which resulted in lower biomass accumulation.

Irrespective of the crop rotation, the average cereal rye biomass produced over the study period followed a strong gradient from north to south (Figure 4). Maximum biomass was produced in the southern areas of the state, and gradually declined northwards. For further analysis, we divided rye biomass data into three distinct groups based on crop rotation to quantify the impact of various weather and soil variables (

Table 3. Quantifying the variability explained (%) by various factors in cereal rye biomass under different crop rotations.

Variables	CRCR (%)	CRSR (%)
		CrC	CrS
Latitude	89 (−773)	67 (−207)	73 (−186)
Precipitation	3 (−38)	1 (−8)	5 (−51)
Temperature	2 (2077)	1 (646)	3 (2257)
Clay	-	7 (33)	5 (118)
Sand	-	2 (19)	-
CEC	-	3 (138)	-
pH	-	1 (701)	-
Residuals	5	17	9

Note: Values in parenthesis represent slope for each variable estimated by developing multiple linear regression model; CRCR = corn-rye-corn-rye; CRSR = corn-rye-soybean-rye; CrC = cereal rye preceding corn; CrS = cereal rye preceding soybean.

). Latitude was found to be the most significant in rye biomass accumulation across the rotations. The significant negative slopes for latitude under all scenarios suggested a decrease in biomass accumulation from south to north, which validated the visual pattern. Our results aligned well with observations from other researchers who found that northern latitudes tend to have a shorter growing season [52], thus, limiting biomass accumulation. Seasonal precipitation (January–April) also expressed a negative slope on cereal rye biomass, such that higher precipitation reduced biomass accumulation. We speculate that this is due to reduced solar radiation received on high precipitation days in the spring, which is common in the U.S Midwest. Soil characteristics, especially clay content, was found to be significant under the CRSR rotation for both phases and recorded a positive influence on rye biomass.

3.3. Soil Organic Carbon

Soil organic carbon simulated by the APSIM model in the top layer (0–45 cm) appeared to have benefitted from the addition of WCCs to the crop rotation. Soil depth up to 45 cm was selected to provide more comprehensive assessment of the SOC quantification and was based on previous studies [53,54]. To quantify the impact of WCCs on carbon sequestration, we analyzed SOC data at 5 year intervals across the whole study area. We estimated the proportion of pixels that fell into different categories based on the observed median values for probability of an increase in SOC (pSOC) and expected value of change in SOC (ESOC) in 2010, 2015, and 2020. pSOC values less than or equal to 0.75 were classified as low representing a low probability of an increase in SOC, whereas pSOC values greater than 0.75 were classified as high probability areas, as suggested by [55]. Similarly, for the ESOC, values lower than or equal to 15% were classified as low expected values, and values greater than 15% were classified as high expected values, based on the difference estimated between WCC and control treatments for SOC by [56]. For the CRCR rotation, 36% of the cropland area in 2010 was found to be under Class 1 when the total land area under this class was reduced to 16% in 2015, and to 9% in 2020 (Figure 5, Table S1). Similarly, the percentage of area under Class 2 followed a similar trend. The area in Class 2 was reduced from 29% in 2010 to 18% in 2015 and was estimated to be around 11% in 2020. The estimated area under Class 3 remained consistently low compared with other classes, starting from 3% in 2010 to 1% in 2020. On the other hand, the area under Class 4 increased from 32% in 2010 to 64% in 2015. In 2020, the area for Class 4 under CRCR rotation was estimated at around 79%, suggesting that at least 79% of croplands in the state of IL have the potential to significantly build up SOC with high probability of success following 16 years of incorporating cereal rye in corn–corn crop rotation. Interestingly, the majority of the area that benefitted from WCC integration under CRCR rotation the fastest was in the southern and western regions of the state which are often associated with lower quality soils.

For the CRSR rotation, 41% of the study area was estimated under Class 1 and 21% under Class 4 in 2010 (Figure 6, Table S2). However, the area under Class 1 reduced to 18% in 2015, and to 7% in 2020, whereas the area under Class 4 increased to 58% in 2015 and to 82% in 2020. The simulated data signify the importance of crop biomass, both WCCs and cash crops, for increasing SOC. The combined area under Class 2 and Class 3 decreased from 38% in 2010 to 11% in 2020. When observing the area under different classes for CRSR, the area that appears to have benefitted most from CRSR is more uniformly distributed from north to south and slightly greater than that seen with CRCR.

On the temporal scale, a large proportion of the area that was converted to Class 4 by the year 2015 was in the western and southern regions of the state. This suggests that WCC benefits with respect to SOC can be harnessed earliest in these regions. We speculate that soil properties, weather conditions, rye biomass, and/or interactions of these variables were largely responsible for such SOC enhancement. By the year 2020, more areas in the central region were classified as Class 4. Most of this new Class 4 area in the central region was previously classified as Class 2. Thus, our results indicate that areas with higher expected change in SOC grow into areas with high probability areas as well, over time. The mean carbon sequestration rate was estimated to be 0.32 Mg C ha⁻¹ yr⁻¹ (1.19 Mg CO₂-eq ha⁻¹ yr⁻¹) for the CRCR rotation and 0.23 Mg C ha⁻¹ yr⁻¹ (0.85 Mg CO₂-eq ha⁻¹ yr⁻¹) for the CRSR rotation, which was similar to those reported by [19].

On further investigation, it appears that for the northern region of the state, CRSR is more beneficial than CRCR with respect to SOC enhancement. Although soybean did not contribute much residue biomass to the soil, it provided a greater growing window for cereal rye to grow in the spring. The CRCR rotation did not allow for such a large growing window for cereal rye in the north, even though corn contributed significantly greater residue biomass compared with soybean. On the other hand, for the cereal rye grown prior to corn, especially in the CRCR rotation, we observed a strong biomass gradient towards the south (Figure 4). Therefore, the net amount of organic matter added to the soil, which included greater rye biomass along with the corn biomass, was much greater than that added during CRSR. Therefore, we speculate that a CRCR rotation offers greater chance of enhancing SOC than CRSR in the southern region of the state.

3.4. Crop Performance

3.4.1. Corn

Unlike SOC, crop performance is more likely to be independent from year to year. Therefore, in our analysis of crop performance, we summarized information across years. For corn, we aggregated the probability of yield increase (pY) and expected yield increase (EY) across all 16 years for the CRCR treatment and across the eight corn years for the CRSR treatment. The aggregated results are presented in Figure 7. Aggregated pY varied from 0.36 to 0.89, and EY varied from −2.4 to 69.7%. Based on the summary of aggregated data, pY less than or equal to 0.5 was categorized as low, and pY greater than 0.5 was categorized as high. Similarly, aggregated EY less than or equal to 5% was classified as low, while EY greater than 5% was classified as high. On average, only 6% of the area fell into Class 1 (low pY and low EY) for corn production when grown under CRCR rotation while Class 4 (dark green) occupied 63% of the area. Classes 2 and 3 on average occupied a combined area of around 31%. For the CRSR rotation, the Class 1 area was estimated to be around 3%, and Class 4 covered around 60% of the area. A combined area of around 37% was classified into Class 2 or Class 3. The results suggest that, irrespective of the crop rotation, on average WCCs benefitted corn production on more than half of the area across the region (high pY and high EY). Similar to SOC, the area that benefitted from the WCC integration was concentrated in the southern part of the state.

Hence, to quantify the impact of various factors on the mean difference in corn yield, we partitioned the variation into various factors using a similar approach as discussed previously (Figure S5). Likewise, we also developed a linear model to assess the magnitude and direction of impact of these factors on mean difference in corn yield between the control group (CC and CS) and treatments (CRCR and CRSR) (Table S3). The findings of variance partitioning suggested that latitude accounts for around 50% of the variability, followed by sand content of the soil, and rye biomass (3% each). Seasonal precipitation (April–October) was found to be a significant predictor. However, it explained less than 1% of the variability. Most of the variables had a positive relationship with the mean difference, except for latitude and seasonal precipitation. The negative slope for latitude suggested that the mean difference for corn yield reduced as we moved from south to north. Similarly, the negative slope of seasonal precipitation implied that the mean difference in yield was alleviated with the increase in seasonal precipitation and was able to compensate for improved soil moisture retention under cover crops.

The YSI provided information about variation in crop yields (yield stability) when WCCs were included in the crop rotations. Around 98% of the area reported YSI greater than 1 for the CRCR rotation when compared with CC rotation reflecting higher yield stability in WCC simulations (Figure 8). Similarly, around 94% of the area recorded YSI greater than 1 for CRSR rotation when compared with CS rotation. The aggregated YSI_CRCR was reportedly higher with a mean of 1.4 ranging from 0.3 to 2.3, compared with YSI_CRSR for eight corn years with a mean of 1.2 ranging from 0.4 to 2.3. The results of YSI analysis suggest that WCCs provide greater stability to corn yield over the long term compared with no cover crops. Contrary to the spatial pattern observed in rye biomass from north to south and SOC, no such trend was observed for yield stability in corn. Hence, when comparing yield stability in WCC rotations versus non cover crop rotations, we speculate that yield stability is not a function of rye biomass.

3.4.2. Soybean

Similar to the corn yield analysis, for classification of the probability of yield increase in soybean (pY), values less than or equal to 0.5 were classified as low probability, whereas values greater than 0.5 were classified as high probability. The average pY across regions varied from 0.44 to 0.59. Similarly, the expected value of yield increase (EY) less than or equal to 5% was categorized as low, and EY greater than 5% as high. The average EY varied from −2% to 15%. Soybean yield benefitted from the inclusion of WCCs as well, even when the WCCs were terminated just one day before planting. When averaged across 8 years, 6% of the soybean production area was categorized as Class 1 and 38% was categorized as Class 4 (Figure S6). The remaining 56% was classified as Class 3. Similar to pY and EY, we did not observe any spatial pattern for YSI in the case of soybean. It is evident from Figure S7, that there is great variation in YSI for soybean when WCCs are grown prior to soybean. The YSI value varied from 0.3 to 1.7 with an average value of 0.87.

Although WCCs terminated just one day prior to soybean did not exhibit any penalties on crop yield, there was no clear pattern observed in the soybean performance due to cereal rye. Part of the reason for an almost random trend could be that we did not consider different termination dates for cereal rye preceding soybean. It may be worthwhile to explore different termination dates for cereal rye grown prior to soybean to better understand its effect on crop performance. On a positive note, a large portion of the state was categorized as Class 3, signifying a higher probability of change in crop yield, even if the expected change was lower. No area was classified under Class 2 (i.e., low probability but high expected change) in the entire study area.

4. Discussion

The resulting posterior densities for the crop parameters optimized for simulations were similar to values reported from prior numerical optimization schemes in the literature [1,37,38]. Although the parameter ‘tt_end_of_juvenile’, which characterizes the thermal time from the end of juvenile stage to terminal spikelet stage, explained significant variability in phenology during GSA, it exhibited substantial variation even after optimization. We speculate this was due to the difference in genotype (G) × environment (E) interaction, and/or the inability of APSIM-Wheat to capture certain underlying processes in cereal rye phenology. Further investigation into unpacking the G × E interactions would require a significantly larger dataset that includes explicit genetic information about different cultivars over diverse environments. In previously conducted studies using APSIM [1,38], a relative root mean square error (RRMSE) of around 56% was reported after calibration when optimizing cereal rye biomass from 2002 to 2014 at a single site in Iowa, U.S. With increased available data, RRMSE of 23% for fall and 56% for spring biomass was reported for the same site in Iowa [38]. In a similar attempt, RRMSE of just 11% for spring biomass was estimated when simulating cereal rye biomass for different planting dates at a single site in Nebraska, U.S. [57]. However, all the above-mentioned studies employed numerical optimization approaches and were limited in terms of number of sites, data types or uncertainty quantification, potentially increasing the risk of overfitting and biased estimation of WCC impacts. Since we applied Bayesian optimization for constraining our model parameters, we were able to generate posterior distributions for crop growth parameters allowing us to quantify the uncertainty around model parameters originating from various cereal rye genotypes and G × E interactions.

Due to the large spatial extent of the simulations, we observed substantial variation, both spatial and temporal, in the cereal rye biomass. However, when cereal rye preceded soybean, significantly greater rye biomass was produced. The underlying process for such variation could potentially be management decisions with respect to cereal rye as discussed by [58], who studied the integration of cereal rye into a corn–soybean rotation in Iowa, U.S. and observed that low temperatures in the fall hindered cereal rye growth. Moreover, most of the biomass accumulation took place in the spring. For the above-mentioned study, soybean was planted within 7 days of cereal rye termination, whereas corn was planted at least 7 days after cereal rye termination. Thus, cereal rye preceding corn on average had two fewer weeks to grow during spring compared with cereal rye preceding soybean.

Additionally, there are several factors that influence SOC when WCCs are integrated into crop rotations, but increased SOC as a function of increased carbon input to the soil (above-ground biomass) is well documented [8]. Growing cereal rye in different combinations for a 2 y corn silage–soybean rotation in Iowa, U.S., enhanced SOM by 15% at 0–5 cm depth, when rye was grown after both corn and soybean in comparison with the control with no rye treatment [59]. In our study, the southern areas of the state produced greater WCC biomass compared with the north regardless of crop rotation, and intuitively those areas reported higher SOC earlier in the simulations than low C input areas. Although SOC enhancement due to WCCs was reported across the state, the time required for significant soil carbon buildup was a major factor as well [60]. Increased SOC due to WCCs is one of the mechanisms by which soil physical properties are improved [61]. This point, however, highlights an important limitation in the APSIM model structure since all soil physical properties including bulk density, hydraulic conductivity, etc., are not dynamically updated as a function of time and management, but are kept constant during the simulation. Therefore, future model improvements are essential for a better representation of change in soil physical properties due to WCC integration.

Cover crop adoption in Illinois remains low, with only ~286,000 ha adopted for cover crops as reported by the 2017 census [62]. The major concern in WCC adoption is that WCCs reduce farm profitability in the short term [63]. These losses can be usually supplemented by various cost-share programs [64] or the biomass can be used for alternative purposes such as cattle grazing [25]. We suggest a precision WCC adoption strategy where we can segregate and quantify the ecosystem services based on different agroclimatic districts of the region [65]. Our research suggests that most of the cropping area in the state of Illinois is responsive to the benefits that come with WCC adoption with potentially negligible yield penalties. However, our research showed that the southern counties in the state tended to show promising results and faster carbon buildup, and may be given priority for cover crop adoption programs. However, results may vary since WCCs require a long term commitment by producers, policymakers, and all other stakeholders. The United Nations Food and Agricultural Organization supported the Lima Paris Action Agenda which aims at increasing the SOC stocks over current stocks by 0.4% annually to achieve sustainable development goals [66]. We did observe a 0.4% annual increase in SOC stocks for our entire study period and the estimated soil carbon sequestration rates (0.15–0.22 Mg C ha⁻¹ yr⁻¹) were higher than those reported by [67]. This indicates that sole adoption of WCCs and no tillage management has the potential to achieve sustainable development goals. However, from a practicality perspective, there may still exist constraints such as sink saturation, or non-permanence of the benefits once the practice is altered [21].

The probability of observing corn yield differences varied significantly across region, as discussed previously. However, we observed an overall positive effect of long term WCC adoption on corn production across IL despite variability in results. When aggregated across years, around 95% of the study area showed positive change in corn yield, when CRCR was compared with CC. Similarly, when CRSR rotation was compared with the CS rotation, around 97% of the area reported a positive effect of incorporating WCC. Our results differ from those obtained by [28], as we did not observe yield penalties for corn by inclusion of non-legume cover crops into the rotation. However, we speculate that terminating cereal rye less than 14 days prior to corn planting may have different effect. Long term application of WCCs improved soil water retention at field capacity when sand content in the WCC treatment was higher than the control [15]. However, the authors credited increased SOC and improved soil aggregation due to WCCs with the improvement in soil water retention. Our analysis suggests that sand content did have a positive relationship with crop performance. Further analysis would be required to distinguish between the effect of sand, soil water content and SOC on yield differences of corn, which is beyond the scope of this study, since we observed increased SOC due to integration of WCCs (cereal rye) as well. Unlike corn, our results are in agreement with those of [28] for soybean, as we did not observe any yield penalty due to non-legume WCCs. Additionally, the current results suggest the need to further explore management options for cereal rye preceding soybean.

Ref. [68] emphasized the importance of adopting a wide array of management practices including no-tillage, cover crops, integrated nutrient management, and precision agriculture to prevent land degradation, while noting that agroforestry and biochar application can provide substantial soil carbon sequestration as alternatives to cover crops. To provide a perspective on such practices, the simulated impact of biochar application in a corn–corn system in Iowa reported a 4% increase in SOC after 30 years [69]. Similarly, a sequestration rate of 0.3 Mg C ha⁻¹ yr⁻¹ with an integrated nutrient management strategy where the recommended dose of nutrients applied through chemical fertilizers was supplemented with 10 Mg ha⁻¹ of farmyard manure, under a soybean-wheat rotation in central India, was reported [70]. The integrated nutrient management demonstrated significantly higher carbon sequestration rate than relying exclusively on chemical fertilizers. Hence, depending solely on WCCs for ecological benefits may not be in best interest of farmers and stakeholders in short-term.

Considering some of the questions that were left unanswered during this study, the focus of future studies should assess the impact of different cereal rye termination dates on crop performance, explore different agroclimatic zones for WCC adoption, and WCC benefits as a function of adoption rate and economic analysis. Apart from the current WCC study, we plan to evaluate the impact of various other cover crops species, such as legumes, using a similar methodology and similar temporal and spatial scales.

5. Conclusions

Our results indicate that a well constrained APSIM model can be used to assess the impact of cereal rye as WCC on a regional level. This work presents the most complete uncertainty accounting of the WCCs’ potential when integrated into common crop rotations (corn-corn and corn-soybean) across the entire cropping area in the state of Illinois. The results of this study indicated that there exists a strong gradient for cereal rye biomass across the state with the north producing lower biomass compared with the south. Additionally, having soybean in the crop rotation provided a longer growing window for the cereal rye, thus, enabling greater biomass accumulation. Irrespective of the crop rotation, the majority of the cropping area was found to be responsive to cereal rye adoption with respect to SOC improvement. However, the SOC enhancement started in the southern parts of the state and followed to the northern parts towards the end of the 16 year simulation period. Corn and soybean yields benefited from WCC integration as well, when compared with non-cover crop rotations. The APSIM modeling suggested that cover crop adoption can provide greater stability to corn production. Cover crops showed great promise in sequestering carbon to tackle climate change and land degradation. Moving forward, more localized research by splitting the entire geographic region into various districts can reveal new insights into WCC adoption at such a large spatiotemporal scale.