Ensemble agroecosystem modeling enhances predictions of crop yields and soil carbon across the United States

Accurately estimating crop yields and soil organic carbon (SOC) dynamics is essential for agricultural planning, carbon accounting, and sustainable land management. However, process-based agroecosystem models often produce divergent estimates due to variations in model structure, parameterization, and underlying assumptions. In this study, we developed a multi-model ensemble framework that integrates three widely used process-based models-Daily Century (DAYCENT), DeNitrification DeComposition (DNDC), and Ecosystem model (ECOSYS)-to simulate crop yields and SOC stock changes (0–30 cm) across cultivated lands of the continental United States (CONUS) at 4 km² spatial resolution. Each model was parameterized using harmonized environmental, soil, and management datasets and evaluated using observed crop yields from the National Agricultural Statistics Service and measured SOC data from the Rapid Carbon Assessment. For the baseline period (2014–2023) under conventional corn–soybean rotation, the ensemble mean showed strong agreement with observations (corn: 7.7 vs. 8.5 Mg ha⁻¹, RMSE = 3.0 Mg ha⁻¹; soybean: 2.5 vs. 3.0 Mg ha⁻¹, RMSE = 1.0 Mg ha⁻¹), while simulated SOC stocks (5.5 vs. 4.8 kg C m⁻², RMSE = 2.5 kg C m⁻²) closely matched measured data. Spatially, the ensemble model projected SOC gains in the Midwest and Southeastern regions and losses in the Great Plains and Western United States, underscoring the importance of region-specific management practices. Overall, the ensemble framework characterizes structural uncertainty by capturing the range of plausible model responses, while the ensemble mean provides a more reliable estimate of system behavior than any individual model. By integrating diverse model perspectives, this approach enhances predictive robustness and reduces dependence on model-specific assumptions. Consequently, it offers a scalable and data-driven framework for assessing soil carbon dynamics and crop productivity across US agroecosystems, while explicitly characterizing and communicating uncertainties associated with model variability.

1 Introduction

Agroecosystems occupy a substantial portion of the global landscape, with agricultural land covering about ∼ 48 million km² (32 %) of the global land area in 2022 (FAO, 2025). These system provide food, feed, fiber, and fuel, and also deliver ecosystem services (ESs) such as soil carbon (C) sequestration, water retention, and nutrient cycling (Robertson et al., 2014). As anthropogenic pressures intensify through consumption patterns, resource demand and weather variability, sustainable management of agriculture practices has become essential to maintain productivity and ecological balance. Enhancing soil organic carbon (SOC) sequestration is one of the key benefits of improved management, contributing to climate solutions and long-term soil health maintenance (Paustian et al., 2016).

Among agricultural practices, crop rotation, particularly the inclusion of legumes such as soybean (Glycine max (L.) Merr.), is widely recognized for its positive effects on SOC accumulation (King and Blesh, 2018; Kumar et al., 2018). Legume-based rotations can improve soil fertility through biological nitrogen fixation (BNF), reduce reliance on synthetic fertilizers, and enhancing carbon inputs from crop residues and root exudates (Drinkwater et al., 1998; Meena and Lal, 2018). They also stimulate microbial activity and stabilize SOC stocks compared to monocultures or less diverse crop rotations (Davis et al., 2012; Xiao et al., 2017). These benefits highlight the importance of understanding how site-specific management strategies and crop rotation influence SOC dynamics across diverse ecoregions and soil conditions.

Among the widely used models are the Daily Century (DAYCENT) model, DeNitrification-DeComposition (DNDC) model, and the Ecosystem Simulation Model (ECOSYS) (Parton et al., 1994; Li et al., 1992; Grant, 1998). Each of these models provides unique structural features and strengths representing plant–soil–atmosphere interactions. DAYCENT a daily time-step version of the CENTURY model simulates C and N fluxes with-in plant–soil systems, incorporating processes such as plant growth, organic matter (OM) decomposition, and soil emissions. It is extensively used for evaluating SOC dynamics and nitrous oxide (N₂O) emissions under different land uses and climate conditions (Del Grosso et al., 2005). DNDC, focuses on coupled C and N biogeochemistry, simulate key emissions including methane (CH₄), carbon dioxide (CO₂), and N₂O (Li et al., 1992). ECOSYS is a more mechanistic model that offers detailed representations of plant physiology, root growth, soil physics, microbial dynamics, and management practices, making it well-suited for capturing complex feedbacks in high-resolution agroecosystem analyses (Grant et al., 2001).

Given the strengths of individual models, there remains substantial variability in predicted crop yields and SOC dynamics due to differences in model structure, parameterization, and underlying assumptions (Tao et al., 2018; Smith et al., 1997). Consequently, recent studies increasingly employ multi-model ensemble approaches to represent uncertainty and enhance predictive reliability. Ensemble modeling leverages the complementary strengths of individual models to characterize a range of plausible system responses under diverse climate and management conditions. This framework provides a more robust basis for environmental assessments by quantifying structural uncertainty while utilizing the ensemble mean to provide a stable estimate of mean agroecosystem conditions (Martre et al., 2015; Basso et al., 2025). The Agricultural Model Intercomparison and Improvement Project (AgMIP) exemplifies this progress, using ensemble to evaluate model sensitivities and parameter uncertainties in crop yield projections (Rosenzweig et al., 2014).

However, most model ensemble studies have focused primarily on crop productivity or provided limited evaluation of SOC dynamics, often constrained by small spatial scales (Rosenzweig et al., 2014; Basso et al., 2025). Furthermore, few studies have combined mechanistic models such as ECOSYS with semi-empirical models like DAYCENT and DNDC, despite their widespread use in soil carbon measurement, reporting, and verification (MRV) and carbon credit calculation. This gap limits comprehensive evaluations of how integrated management strategies influence both agronomic outcomes and ecosystem services (ESs).

In this work, we develop and benchmark an ensemble agroecosystem modeling framework integrating DAYCENT, DNDC and ECOSYS to predict crop yields (corn and soybean) and SOC dynamics across the continental United States (CONUS). This framework bridges agronomic, environmental, and computational disciplines, enabling improved carbon accounting and sustainability assessments at continental scales. By combining models with distinct process representations, we assess both consistency and divergence of predictions and identify areas of model convergence that inform robust management strategies. Unlike previous efforts, this study emphasizes the spatial variability in multi modal agreement across US croplands, providing valuable insights into the design of resilient agroecosystems under weather variability. Furthermore, the results provide a spatially explicit baseline that can inform policymakers, and stakeholders in prioritizing regions for carbon MRV and identifying areas at higher risk for SOC decline under current crop rotations.

The specific objectives of this study were to: (a) develop an ensemble agroecosystem modeling framework for CONUS croplands; (b) simulate crop yields and SOC stock changes under baseline conditions from 2014 to 2023; and (c) compare model outputs from DAYCENT, DNDC, and ECOSYS to characterize variability and uncertainty in predicted yields and SOC stocks. We hypothesized that: (1) the ensemble of DAYCENT, DNDC, and ECOSYS would capture the uncertainty of SOC stock change estimates; and (2) ensemble mean corn and soybean yields would align more closely with observations while exhibiting spatially heterogeneous responses linked to geographic and environmental gradients.

2 Methods

2.1 Study area and model input

This study encompassed three primary land cover types (croplands, pastures, and grasslands) identified using National Land Cover Database (NLCD) (USGS, 2023). The daily precipitation and minimum and maximum air temperature data for the 30-year period (1994–2023) were obtained from Daily Surface Weather and Climatological Summaries (DAYMET) (Thornton et al., 2016). The key soil parameters for model initialization included multilayer soil texture, bulk density, soil moisture retention characteristics, hydraulic conductivity, soil organic carbon (SOC), and pH, derived from the Soil Survey Geographic Database (SSURGO) (Soil Survey Staff, 2023). The models (DAYCENT, DNDC, and ECOSYS) represent vertical soil processes using discretized soil layers. Within each layer, key properties such as soil texture (sand and clay fractions), bulk density, pH, and hydraulic characteristics are specified, enabling simulation of depth-dependent dynamics of soil carbon, water, and nutrient processes. The study simulated a corn–soybean rotation spanning 2014 to 2023 to evaluate crop management impacts on soil biogeochemical processes. Nitrogen (N) fertilizer was applied at a rate of 185 kg N ha⁻¹ during the corn year following planting. Tillage operations included a spring chisel plow of 30 cm and a fall field cultivator to 20 cm, representing conventional management practices. In contrast, soybean cultivation followed a no-till (NT) approach with minimal soil disturbance during the growing season. We selected a 4 km spatial resolution to balance between computational efficiency and the need to adequately capture large-scale spatial variability across the CONUS study domain. This resolution is appropriate for regional to national-scale assessments of crop productivity and soil carbon dynamics, including broad carbon accounting applications. However, it may not fully capture fine-scale heterogeneity in soil properties, management practices, and microclimatic conditions that influence carbon fluxes at the field level. This rotation setup reflects typical Midwestern US agricultural practices and provides a consistent framework for evaluating long-term SOC dynamics and crop productivity across multiple models. The schematic overview of the modeling framework and workflow is presented in Fig. 1.

Figure 1 Schematic representation of the workflow for large-scale application of the ensemble agroecosystem model (DAYCENT, DNDC, ECOSYS). Environmental and management inputs were derived from multiple datasets, including Daymet (daily surface weather data), SSURGO (Soil Survey Geographic Database), NASS (National Agricultural Statistics Service crop yield surveys), RaCA (Rapid Carbon Assessment soil carbon database), and NLCD (National Land Cover Database).

2.2 Process model description and setup

The three models selected for this study, DAYCENT, ECOSYS, and DNDC, were chosen because they represent a broad spectrum of biogeochemical modeling frameworks, ranging from semi-empirical pool-based logic to high-complexity mechanistic kinetics. These models are widely utilized in regional carbon assessments and national greenhouse gas inventories, and their inclusion allows for a robust comparison of how different process representation of soil-plant interactions influence regional SOC and yield projections.

DAYCENT is a daily time step, process-based agroecosystem model developed to simulate the biogeochemical cycles at a point scale (Del Grosso et al., 2002). It is daily version of the CENTURY model (Parton et al., 1998). It can simulate exchange of carbon and soil nutrients (N, phosphorus, P, and sulfur, S) between the atmosphere and terrestrial ecosystems, along with soil water and temperature dynamics (Schimel et al., 2001). In the DAYCENT model, SOC is predicted based on OM decomposition and C cycling among multiple soil pools, including surface litter (structural and metabolic), and fast, slow and passive C pools (Parton et al., 1998; Del Grosso et al., 2002). Carbon inputs arise from plant residues and organic amendments, while decomposition is driven by soil temperature, moisture, and texture. DAYCENT simulates daily C fluxes between soil pools, accounting for mineralization, immobilization, and heterotrophic respiration losses (Parton et al., 1994). Additionally, DAYCENT also represent various management practices including tillage, crop rotation, residue retention, and fertilization. Crop yield is simulated through net primary production (NPP) as a function of sunlight, temperature, soil water, and nutrients. The assimilated carbon in partitioned among plant components including grain, roots and shoots based on their C:N ratios. The final yield is determined using the harvest index, which defines the proportion of aboveground biomass converted into harvestable product, such as grain. This harvest index is crop-specific and sensitive to environmental and stress conditions, including water limitation (drought), nutrient stress (particularly nitrogen deficiency), and temperature extremes (heat or cold stress). This sensitivity enables DAYCENT to represent realistic interactions among biophysical processes and management practices. The model was initialized with a long-term spin-up (1000 to 4000 years) to achieve the steady-conditions for soil nutrient pools, representing pre-disturbance natural land conditions prior to human intervention. For the equilibrium simulations, ecoregion-specific native grass were grown, based on data from the US Environmental Protection Agency Inventory of US GHG emissions and sinks (USEPA, 2015). Historical land use and land cover data were compiled from historical datasets (Ogle et al., 2010). Historical management practices from initial tillage to the modern agricultural period was compiled at the Major Land Resource Area (MLRA) level from multiple historical data sources (Ogle et al., 2010). Contemporary agricultural management was represented based on the multiyear analysis of the National Land Cover Database. Additionally, historical database was used to estimate the average annual nitrogen fertilizer rates for each agricultural crop across different ecoregion (USEPA, 2005).

DNDC is a process based soil biogeochemical model developed to simulate C and N cycling in agricultural ecosystems (Li et al., 1994, 1992). It integrates key soil physical and chemical properties (e.g., soil texture, bulk density and soil pH) to represent soil biogeochemical processes. DNDC also includes a detailed plant growth module that simulates the entire crop life cycle, linking plant development with soil C and N dynamics. The model partitions the net assimilated C among leaves, stems, roots, and grain according to the crop's phenological stage. Final crop yield is determined at the end of the simulation as the total biomass allocated to the grain at physiological maturity. Due to its explicit representation of plant growth and decomposition processes, DNDC is widely used to evaluate the effects of agricultural management practices on crop yields, SOC dynamics and soil emissions (Ding et al., 2023; Jiang et al., 2021).

ECOSYS is a process-based ecosystem model grounded in biophysical and biochemical principles that integrates water, energy, C and nutrient cycles (Grant et al., 2001). It explicitly represents interactions among soil, plant, and atmospheric components with high mechanistic details. Due to its comprehensive process representation, including plant photosynthesis and hydraulics, dynamic C and nutrient allocations in plants and microbes, soil biogeochemistry, and agricultural management practices, ECOSYS is a robust tool for simulating agroecosystem responses to diverse environmental and management conditions (Zhou et al., 2021; Grant et al., 2006; Shekoofa et al., 2021). The model incorporates detailed physical, chemical and biological processes to simulate plant-microbe-soil-atmosphere interactions across diverse soil types, climatic regimes, and land-use practices. It explicitly models nutrient exchange between microbes and plants, allowing nutrient-based evaluation of management effects on crop productivity and SOC dynamics. Crop yield is calculated based on dynamic plant C allocation, which balances photosynthetic C gain and respiratory losses, and partitions assimilated biomass. Grain filling is regulated by canopy temperature, available carbon and nutrient reserves, and physiological constraints, allowing ECOSYS to capture crop yield responses to environmental stress and management variability. The models handle vertical soil discretization according to their distinct structural designs, aligned with our primary evaluation depth of 0–30 cm.

2.3 Model Evaluation

Each of the three models was calibrated using county level crop yield data from the National Agricultural Statistics Service (NASS) (NASS, 2025) and site-level SOC data from the Rapid Carbon Assessment (RaCA) dataset (Wills et al., 2014) across CONUS. Model parameters were optimized to minimize the root mean square error (RMSE) and normalized RMSE (% RMSE). Relative metrics like the Nash–Sutcliffe Efficiency (NSE) and R² were intentionally excluded, as they can be misleadingly inflated in large-scale regional contexts where model predictions closely align with a strong regional observation mean (Gupta et al., 2009; Klotz et al., 2024).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(P_i-O_i{)}^{\mathrm{2}}}{n}}}\text{(2)}& {\displaystyle }& {\displaystyle }\text{NRMSE}={\displaystyle \frac{\sqrt{\frac{{\sum }_{i=\mathrm{1}}^n(P_i-O_i{)}^{\mathrm{2}}}{n}}}{\widehat{O}}}\times \mathrm{100}\end{array}$$

Figure 2 Comparison of observed crop yields (Mg ha⁻¹; NASS, 2014–2023) with simulated yields and root mean square error (Mg ha⁻¹) from individual agroecosystem models and the ensemble approach for (a) corn and (b) soybean.

To explicitly evaluate the risk of model overfitting and assess framework stability prior to final execution, a stratified data-splitting protocol was implemented using the createDataPartition framework in the caret package. Continuous datasets (RaCA SOC and NASS yields) were partitioned into equal-sized numeric bins based on quantiles, and 50 % of the data points were randomly sampled for calibration, leaving the remaining 50 % as an independent validation subset. This cross-validation protocol yielded identical error distributions between the independent groups specifically, an ensemble SOC RMSE of 4.7 kg m⁻² (calibration) vs. 4.6 kg m⁻² (validation), a corn yield RMSE of 3.0 Mg ha⁻¹ (calibration) vs. 3.0 Mg ha⁻¹ (validation), and a soybean yield RMSE of 1.0 Mg ha⁻¹ (calibration) vs. 1.0 Mg ha⁻¹ (validation). Because this identical performance demonstrated that the modeling framework was structurally stable and robust against overfitting, the data subsets were pooled for the final simulations. This approach maximizes geographic continuity across the CONUS domain and fully utilizes the available high-resolution datasets. This approach provides a robust upper bound of current process-based modeling capabilities at a 4 km² scale, representing the collective central tendency of diverse mechanistic assumptions.

Where, P_i and O_i are the predicted and observed values, respectively, and n is the total number of observations. After model evaluation, each model was updated using a calibrated parameter set, and large-scale simulations were executed on the Sandia National Laboratories high performance computing (HPC) system in a parallel environment.

3 Results

3.1 Model Evaluation of crop yield and SOC

A comparison of simulated corn yields across the CONUS croplands from three agroecosystem models and their ensemble median against observed yields from the NASS dataset showed that the ensemble model provided the closest agreement with observations (Fig. 2a). The ensemble exhibited the lowest variability and achieved the smallest root mean square error (RMSE = 3.0 Mg ha⁻¹, NRMSE = 39 %) (Fig. 2a). Spatial distribution of RMSE for corn ranged from approximately 0–9.4 Mg ha⁻¹ across models (Fig. S1a in the Supplement). Its median corn yield of (7.7 Mg ha⁻¹) closely matched the observed median of 8.5 Mg ha⁻¹. Among individual models, DAYCENT and DNDC performed moderately well, producing RMSEs of 3.8 and 5.0 Mg ha⁻¹, and median yields of 9.4 and 8.7 Mg ha⁻¹, respectively. In contrast, ECOSYS exhibited the greatest deviation, with an RMSE of 4.1 Mg ha⁻¹ and a median yield of only 4.4 Mg ha⁻¹, reflecting a consistent tendency to underpredict yields and higher variability across the CONUS region (Fig. 2a).

For soybean yields, the ensemble model again demonstrated the best performance, yielding the lowest root mean square error (RMSE = 1.0 Mg ha⁻¹, NRMSE = 41 %) and most accurate representation of the observed yield distribution (Fig. 2b). The ensemble median soybean yield (2.5 Mg ha⁻¹) was in close agreement with the observed median (3.0 Mg ha⁻¹). Among the individual models, DAYCENT and DNDC yielded RMSEs of 1.3 and 2.1 Mg ha⁻¹ with median yields of 3.2 and 2.7 Mg ha⁻¹, respectively (Fig. 2b). ECOSYS again underpredicted soybean yields, with an RMSE and a median yield of 1.3 Mg ha⁻¹. Overall, the ensemble approach minimized prediction errors and more accurately represented the central tendency of observed yields compared with individual models. These results highlight the advantage of multi-model ensemble frameworks in enhancing prediction accuracy and reducing uncertainty in large-scale crop yield simulations. Spatial distribution of RMSE ranged from approximately 0–2.9 Mg ha⁻¹, indicating closer agreement with NASS observations (Fig. S1b).

A comparison of simulated SOC stocks (0–30 cm) from three agroecosystem models (DAYCENT, ECOSYS, and DNDC) and their ensemble against observations from the RaCA dataset revealed distinct differences in model performance (Fig. 3). The RaCA observations indicated a median SOC stock of 4.8 kg C m⁻² across CONUS for 0–30 cm layer. Among the individual models, DAYCENT performed most consistently with observations, producing a median SOC of 5.5 kg C m⁻² and an RMSE of 2.7 kg C m⁻². ECOSYS underestimated SOC, yielding a median of 3.2 kg C m⁻² and an RMSE of 3.0 kg C m⁻², suggesting limitations in its representation of SOC accumulation processes. In contrast, DNDC consistently overestimated SOC, with median of 7.7 kg C m⁻² and the largest RMSE (4.9 kg C m⁻²) among the three models (Fig. 3). The ensemble approach effectively reduced model-specific biases, producing a median SOC estimate of 5.5 kg C m⁻² with the lowest RMSE (2.5 kg C m⁻²) across CONUS (Fig. 3). Collectively, these results demonstrate that ensemble modeling enhances agreement with observed data and improves the reliability of SOC stock estimates at regional and national scales. Spatially, the ensemble provides more consistent predictions in regions where individual models diverge, highlighting its strength in capturing heterogeneous environmental and management conditions.

Figure 3 Comparison of observed soil organic carbon (SOC) stocks (0–30 cm) (kg C m⁻²) from the RaCA dataset with SOC simulated and root mean square error (kg C m⁻²) by individual agroecosystem models and the ensemble approach.

Figure 4 Spatial distribution of annual corn yield (Mg ha⁻¹) across the CONUS simulated by the ensemble and individual agroecosystem models (DAYCENT, DNDC, and ECOSYS) for the 2014–2023 period.

3.2 Model prediction of corn yield

During the 10-year simulation period (2014–2023), DAYCENT estimated annual corn yields ranging from 1.0 to 15.3 Mg ha⁻¹, with a mean of 8.8 Mg ha⁻¹. DNDC projected yields between 1 and 15 Mg ha⁻¹ (mean = 8.7 Mg ha⁻¹), whereas ECOSYS computed a similar range but a substantially lower mean of 5.95 Mg ha⁻¹ (Fig. 4). All three models captured the broad spatial pattern of high corn productivity across the US Midwest Corn Belt, particularly in central Iowa, Illinois, and Indiana. Regional difference among models reflected contrasts in structure representation and parameterizations. DNDC projected relatively higher yields in the Mississippi River Basin and the Southeast US, while DACYENT simulated more moderate yields in the same regions (Fig. 4). By contrast, ECOSYS estimated comparatively higher productivity across parts of the Northern Plains and Northwestern US The ensemble yield map effectively moderated these model-specific extremes, generating a spatially coherent pattern that preserved agreement within high-yield zones while diminishing the influence of outlier projections. Overall, the ensemble framework produced robust regional yield estimates by integrating model variability, leveraging complementary strengths, and redacting individual model biases.

Figure 5 Spatial distribution of annual soybean yield (Mg ha⁻¹) across the CONUS simulated by the ensemble and individual agroecosystem models (DAYCENT, DNDC, and ECOSYS) for the 2014–2023 period.

3.3 Model prediction of soybean yield

From 2014 to 2023, DAYCENT simulated annual soybean yields ranging from 0.1 to 5 Mg ha⁻¹ with a mean of 2.8 Mg ha⁻¹. DNDC projected yields within a similar range (0.1–5 Mg ha⁻¹) and a mean of 2.7 Mg ha⁻¹, while ECOSYS estimated a slightly broader range (0.01 to 6 Mg ha⁻¹) with a mean of 2.5 Mg ha⁻¹ (Fig. 5). All three models reproduced the general spatial pattern of highest soybean productivity across the Midwestern and central US, particularly in major soybean producing regions such as Iowa, Illinois, and Indiana. Model-specific variability emerged in both the spatial extent and magnitude of yield predictions. DNDC projected high yields concentrated in central Iowa and surrounding areas, whereas DAYCENT and ECOSYS predicted more spatially uniform yield distributions across the Corn Belt (Fig. 5). Additionally, DAYCENT simulated comparatively higher yields across parts of the Southeast US.

The ensemble yield map integrated outputs from all three models, producing a harmonized spatial distribution characterized by a concentrated high-yield zone in the central Midwest and a gradual decline in productivity toward coastal and arid western regions. This ensemble approach provided a more coherent and robust estimate of soybean productivity across diverse agroclimatic zones, effectively smoothing local extremes and reducing model-specific outliers.

Figure 6 Spatial distribution of simulated soil organic carbon (SOC) stock changes ($\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) over a 10-year period (2014–2023) using individual agroecosystem models and the ensemble approach.

3.4 Model prediction of soil organic carbon change

Figure 6 illustrates the spatial distribution of projected annual SOC change ($\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) under a 10-year (2014–2023) conventional-till corn–soybean rotation across CONUS, based on simulations from the three agroecosystem models and their ensemble mean. Positive values represent SOC sequestration, whereas negative values indicate net C emissions. DAYCENT consistently projected widespread SOC sequestration across the Midwest and eastern US, with rates often exceeding 0.02 $\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in many areas (Fig. 6). In contrast, DNDC and ECOSYS exhibited greater spatial heterogeneity, including notable zones of SOC decline in the Northern Great Plains and parts of the Midwest. DNDC simulated SOC losses surpassing 0.01 $\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, particularly within high-yielding zones (Fig. 6). The ensemble map integrated these model-specific responses, capturing the broad sequestration trends indicated by DAYCENT while incorporating SOC-loss zones highlighted by DNDC and ECOSYS. Figure S2 in the Supplement illustrates the spatial distribution of RMSE between ensemble model-simulated SOC and RACA observations across the conterminous United States, highlighting regional variability in model performance. Overall, the ensemble projection suggested a modest net increase in SOC across much of the Corn Belt, underscoring the capacity of conventional corn–soybean systems to enhance carbon sequestration under region-specific environmental and management conditions. However, these estimates represent modeled outcomes and should be interpreted with caution, as they may vary depending on site-specific conditions and management practices.

Table 1 Summary of biogeochemical process representations and model structures for the three-model ensemble.

Note: Daycent = Daily Century; DNDC = Denitrification-Decomposition; Ecosys = Ecosystem model.

4 Discussion

The ensemble modeling approach produced corn yield estimates across the CONUS that closely aligned with the observed data, underscoring the value of multi-model ensembles in reducing individual model biases and improving prediction accuracy (Asseng et al., 2013; Martre et al., 2015; Basso et al., 2025). The ensemble model benefited from complementary model structures, error compensation and statistical averaging, resulting in enhanced stability and performance across diverse environmental gradients. The moderate accuracy of DAYCENT and DNDC was consistent with prior evaluations demonstrating their strengths in simulating corn productivity under temperate conditions, particularly in the US Midwest (Del Grosso et al., 2002; Campbell et al., 2014; Ingraham and Salas, 2019). In contrast, ECOSYS substantially underpredicted corn yield, likely due to its complex representation of coupled root and canopy processes, which increases sensitivity to input uncertainty and parameterization limitations in large-scale applications (Liu et al., 2024). In our large-scale application, the use of standardized N-input data likely triggered disproportionate N-stress within the model, leading to lower-than-observed yields. ECOSYS employs a detailed rhizosphere sub-model where N-uptake is governed by ion diffusion and mass flow this can lead to localized N-stress that is more extreme than observation (Grant et al., 2001). Regional anomalies, such as DNDC's higher corn yield projections in the Southeast, are likely influenced by its representation of water stress and nutrient availability under less favorable climatic conditions, as well as its relatively linear response of crop yield to precipitation (Kaur et al., 2023). Collectively, these results supported the hypothesis that ensemble predictions provide more reliable and robust estimation of corn yield across diverse agroecological zones in CONUS.

Similarly, the ensemble model demonstrated superior performance in predicting soybean yield, achieving the lowest RMSE (1.0 Mg ha⁻¹) and reproducing the observed median yield (2.5 Mg ha⁻¹). This finding aligns with previous research highlighting the advantages of ensemble modeling in stabilizing projections across variable cropping systems (Battisti et al., 2017). The relatively strong performance of DAYCENT and DNDC reflected their ability to simulate soybean physiology and biological nitrogen fixation (BNF) processes in temperate regions. However, DNDC's spatial overestimation in Iowa is likely resulted from its limited representation of soil drainage and water balance (Helmers et al., 2012). ECOSYS underpredicted soybean yields, possibly due to its sensitivity to soil moisture dynamics and restricted calibration in drier regions. These model-specific biases highlight the ongoing challenges in accurately representing complex interactions among soil texture, weather variability, and cultivar adaptability, further reinforcing the value of ensemble approaches for improving yield estimation.

Ensemble projections of SOC stock changes (0–30 cm) across CONUS under a 10-year corn–soybean rotation indicate a modest net sequestration. In the Corn Belt, SOC changes exhibit substantial spatial variability, with both gains and losses (approximately −2 to +2 $\mathrm{t}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), resulting in a relatively small positive average across the region. This pattern is broadly consistent with previous studies that report SOC sequestration potential under specific management or rotation scenarios (Grace and Robertson, 2021; Wu et al., 2025), although those studies often consider alternative practices or crop shifts not explicitly represented here. The magnitude of the modeled SOC increase should be interpreted in the context of model uncertainty, as reflected by the ensemble's tendency to overestimate SOC relative to the RaCA dataset. Therefore, while the results suggest a potential for modest carbon gains, their significance at regional scale remains sensitive to both spatial heterogeneity and model accuracy. Strong performance of DAYCENT in capturing SOC accumulation within the Midwest aligns with its established capacity to simulate long-term C cycling and residue turnover in US cropping systems (Basso et al., 2025). In contrast, DNDC and ECOSYS underpredicted SOC stocks, likely due to model sensitivities to decomposition rates, soil moisture retention, and microbial activity, which are critical determinants of SOC turnover. The greater variability observed with DNDC was consistent with prior reports attributing its coarser parameterization of SOC pools and limited regional calibration, as its original design emphasis on CH₄ and N₂O emissions rather than SOC dynamics (Zhang et al., 2014). The ensemble's mean SOC and reduced RMSE illustrate its effectiveness in mitigating structural uncertainty and reproducing the SOC patterns captured in the RaCA database. These results underscore the promise of ensemble-based approaches for improving national scale carbon accounting and guiding sustainable land management strategies.

The divergent SOC trajectories across the models shows how different biogeochemical model interpret the management impacts differently. The widespread SOC decrease projected by DNDC, particularly in high-yielding zones, is driven by its Double Monod kinetic framework, which simulates microbial growth as a function of simultaneous carbon and mineral nitrogen availability (Li et al., 1992). In this framework, abundant nitrogen leads to mineralization of existing organic matter to satisfy stoichiometric demands faster than humification of crop residue. Conversely, DAYCENT employs a conservative pool-based approach where nitrogen primarily boosts plant productivity and residue inputs rather than accelerating the decay of stabilized pools (Parton et al., 1998). ECOSYS uses microbial carbon use efficiency and physical protection within soil aggregates, assuming that high-input carbon is stabilized as microbial biomass rather than respired (Grant, 1998).

This study uses a uniform set of crop phenotypic and parameter values for corn and soybean across the entire spatial domain. While this standardization enabled large-scale simulations and computational efficiency, it inherently omitted genotype–environment interactions and region-specific cultivar adaptations. In ECOSYS, for example, crop phenology is highly sensitive to thermal time and maturity group parameters, both of which vary significantly across US regions. Consequently, the underperformance of ECOSYS relative to DAYCENT and DNDC likely reflects a sensitivity penalty where the model's mechanistic complexity was constrained using simplified, uniform maturity groups. Previous work have demonstrated that regional calibration of maturity group is essential for accurately reproducing phenological development (Li et al., 2022); thus, using a single maturity group may have misrepresented region-specific growth dynamics and yield responses. Despite this, ECOSYS was retained in the ensemble to ensure the structural uncertainty range accounts for high-complexity mechanistic feedback such as plant–water relations and microbial kinetics that simpler, pool-based models might overlook. Although this approach supports continental scale consistency, it introduces limitations in representing localized cultivar adaptation and environmental responses, highlighting an important direction for future ensemble model refinement. Additionally, the 4 km spatial resolution, while appropriate for regional-scale assessments, may not capture fine-scale variability in soils and management practices, which could limit its applicability for field-scale applications such as carbon credit quantification. Finally, the limited ensemble size in this study may not capture the full structural uncertainty found in larger initiatives like AgMIP, as additional models could introduce different sensitivities that would likely widen the projected uncertainty bounds.

5 Conclusions

This study advances multi-model integration in agroecosystem science by developing and benchmarking an ensemble framework that combines three process-based models, DAYCENT, DNDC, and ECOSYS, to improve large scale predictions of crop yields and soil carbon dynamics across the continental United States. The ensemble approach outperformed individual models, producing lower root mean square errors (RMSEs) and stronger agreement with observed corn and soybean yields and SOC stocks in terms of both mean values and variability. Ensemble results captured the central tendency of observed data more effectively, while generating spatially coherent yield patterns. For SOC dynamics, ensemble projections represented both sequestration potential and emission hotspots, capturing net carbon gains in the Midwest while reducing over and under estimation patterns from DNDC and ECOSYS, respectively. Overall, this study highlights the value of an ensemble modeling framework for enhancing regional to national-scale assessments of crop productivity and soil carbon dynamics. By characterizing structural uncertainty and utilizing the ensemble mean to mitigate individual model biases, this approach provides a more robust scientific basis for carbon measurement, reporting, and verification. However, its application to farm-level carbon credit quantification and specific regenerative agriculture initiatives should be considered preliminary, as these applications may require further model refinement and higher-resolution spatial data to capture local-scale variability. Ensemble agroecosystem modeling provides a critical foundation for advancing evidence-based land management and sustainable agricultural policy in a changing environmental condition.