A food crop yield emulator for integration in the compact Earth system model OSCAR (OSCAR-crop v1.0)

This paper presents the development, validation, and preliminary application of a sub-national scale crop yield emulator to be integrated into the compact Earth system model OSCAR. The emulator simulates yields for four major food crops: maize, rice (two growing seasons), soybean, and wheat (spring and winter varieties), in alignment with the Agricultural Model Intercomparison and Improvement Project (AgMIP) and the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) framework. Key drivers include atmospheric CO₂ concentration (represented as C), growing season temperature (T), water availability (W), and nitrogen fertilization (N). The emulator is trained on an ensemble of process-based crop model simulations from AgMIP's Global Gridded Crop Model Intercomparison Projects (GGCMI), which is based on the ISIMIP Phase 3 protocol. These crop models used bias-corrected historical and future (SSP126, SSP370, and SSP585) climate scenarios under fixed human direct forcing to estimate yield responses to C, T and W. Evaluation of the emulator against the crop model outputs demonstrates the emulator's ability to replicate complex model behavior with high fidelity. Additionally, the emulator-derived yield sensitivities to CO₂ and temperature are consistent with those observed in field experiments, reinforcing its empirical robustness. Historical simulations incorporating time-varying nitrogen inputs show significantly improved agreement with FAO yield statistics, underscoring the emulator's reliability over the historical period and its potential for future impact assessments. This study provides a computationally efficient yet empirically grounded tool for representing crop yield responses, providing a middle-ground between complex crop models and statistic models. The developed crop emulator facilitates probabilistic projections across large ensembles of climatic and socio-economic scenarios at policy-relevant, sub-national scales. Potential applications include integrated assessments of future food security under climate and land-use change, as well as evaluations of bioenergy with carbon capture and storage (BECCS) potential from crop residues.

1 Introduction

As the global population grows and climate change poses increasing risks to agricultural productivity, ensuring food security becomes increasingly critical and challenging. Robust projections of major food crop yields under future climate conditions are essential. Crop models are fundamental tools for such projections and are broadly classified into statistical models and process-based models. Statistical models (Lobell and Asseng, 2017) largely built upon historical relationships between climatic variables and crop yields. However, their predictive power diminishes under future high-emission scenarios where conditions are outside of the historical range (Ciscar et al., 2018). In contrast, process-based models simulate the biophysical processes underlying crop growth and responses to management and climate factors (Franke et al., 2020a; Di Paola et al., 2016). However, their computational cost limits their application to large scenario ensembles. To address this trade-off between complexity and computational efficiency, crop yield emulators have emerged as a promising alternative (Abramoff et al., 2023; Blanc, 2017; Folberth et al., 2025; Franke et al., 2020a; Liu et al., 2023; Ringeval et al., 2021). These emulators approximate the behavior of more complex crop models, thereby enabling large-ensemble simulations for robust impact assessments (Folberth et al., 2025).

The Global Gridded Crop Model Intercomparison Project (GGCMI), part of the Agricultural Model Intercomparison and Improvement Project (AgMIP) (Rosenzweig et al., 2013), provides harmonized projections of food crops across diverse climate scenarios (Franke et al., 2020a; Jägermeyr et al., 2021; Müller et al., 2017). GGCMI Phase 1 enabled the development of a phenology-based crop emulator for potential yields (Ringeval et al., 2021), and GGCMI Phase 2 facilitated the development of a crop emulator trained on historical weather inputs with perturbations to CTWN at grid cell level (Franke et al., 2020a; Müller et al., 2024). In GGCMI Phase 3, aligned with the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (Frieler et al., 2024), crop yield simulations are generated using the latest generation of complex crop models forced by bias-adjusted CMIP6 climate simulations (Jägermeyr et al., 2021). Despite the wide range of outcomes due to different model structures, parameterization schemes, calibration processes and input data quality (Folberth et al., 2019; Müller et al., 2024), the projections in GGCMI Phase 3 exhibit reduced uncertainty for rice and soybean and enhanced robustness for maize and wheat (Jägermeyr et al., 2021). Therefore, a sub-national crop emulator calibrated against the state-of-the-art GGCMs and integrating information from long-term field experiments is developed in this study. Through a multi-ensemble approach, the crop emulator covers inter-model uncertainties from both climate and crop models and provides probabilistic distributions of crop yields rather than deterministic yields, improving the reliability of future yield projection (Maiorano et al., 2017).

This paper provides a comprehensive description of this crop emulator, named OSCAR-crop v1.0, regarding its design, calibration, and validation. Section 2 outlines the emulator framework and data preprocessing procedures. Section 3 details the calibration methodology and driver-specific yield responses. Section 4 evaluates model performance through both in-sample (ISIMIP3b) and out-of-sample (ISIMIP3a) validation against original GGCMs' results. Section 5 presents a comparative analysis of emulator-derived sensitivities against experimental data, and the comparison of simulated crop yields with FAO national records to further validate the emulator. Finally, Sect. 6 synthesizes key findings, discusses uncertainties, and presents potential future applications.

2 Model overview and preprocessing

2.1 Overview

The OSCAR-crop v1.0 is developed at a regional level with annual temporal resolution to facilitate region-specific applications. The emulator encompasses 311 regions in total, providing sub-national modelling capabilities for six large-area countries – Australia, Brazil, Canada, China, Russia, and the USA – to capture internal spatial heterogeneity. For the remainder of the globe, the model operates at a national level. The emulator incorporates maize (mai), rice from the first and second growing seasons (ri1 and ri2), soybean (soy), and two types of wheat (spring wheat and winter wheat, swh and wwh). Its design aligns with the probabilistic framework of the compact Earth system model OSCAR, enabling robust numerical projections for large ensembles of climatic and socio-economic scenarios (Gasser et al., 2017; Quilcaille et al., 2023). The model offers two operational modes: a standalone application requiring external CTWN inputs, and integration within a simple climate model (SCM) framework. For integration with OSCAR, the crop emulator acts as a module, receiving its CTW inputs from the host. This coupling requires the SCM to provide necessary variables to drive the emulator. Moreover, it allows for easy adaptation to coarser regional scales, making it a valuable tool for integrated assessment modelling (Ruane et al., 2017).

Crop yields are influenced by a range of factors, including temperature (Wang et al., 2017), water availability (Proctor et al., 2022), solar radiation (Laub et al., 2022), ambient CO₂ (Toreti et al., 2020), soil quality (Qiao et al., 2022), and human management (Ahvo et al., 2023). Four key drivers including atmospheric CO₂ concentration ([CO₂]), growing season temperature, water input, and nitrogen fertilization are incorporated to streamline the crop emulator, representing critical climatic and management variables. All crops are classified into two categories based on the irrigation system: rainfed (noirr) and fully irrigated (firr). The core outputs of the emulator are crop yields. The structure of the emulator is shown in Fig. 1.

Figure 1 Crop yield emulator structure.

For illustration purposes only, a tiered aggregation strategy is employed to improve the presentation of regional results in this study, drawing attention to the significance of major-producing regions. First, sub-national regions are aggregated to the national level, and EU27 countries are combined into a single entity. Subsequently, the five leading crop-producing regions (R1 to R5), identified based on current production data (Fig. S2 and Table S1 in the Supplement), are presented individually for each crop, while all other regions are collectively aggregated (R6).

2.2 Data preprocessing

The emulator development utilizes data from the ISIMIP repository (https://data.isimip.org/search/, ISIMIP repository, 2024). The raw input variables provided in the repository encompass daily temperature and precipitation, annual nitrogen fertilizer inputs, crop-specific land areas, and crop calendars, all at the grid cell level. Raw output variables include crop yield responses and crop yields simulated by eight state-of-the-art GGCMs, including CYGMA1p74 (Iizumi et al., 2017), EPIC-IIASA (Balkovič et al., 2014), ISAM (Gahlot et al., 2020), LDNDC (Haas et al., 2013), LPJmL (von Bloh et al., 2018; Lutz et al., 2019), PEPIC (Liu et al., 2016), PROMET (Hank et al., 2015; Mauser et al., 2015; Zabel et al., 2019), and SIMPLACE-LINTUL5 (Webber et al., 2018). These GGCMs are run under pre-industrial (picontrol), historical, and future (SSP126, SSP370, and SSP585) climate scenarios, projected by five Earth system models (ESMs), namely GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0 and UKESM1-0-LL. The default experiments in ISIMIP3 fix direct human forcing (e.g., land-use, nitrogen fertilizer use) at the 2015 level (Frieler et al., 2024), while [CO₂] and climate variables vary according to the specified climate scenario. To specifically assess the impact of CO₂ on long-term crop yield projections, sensitivity experiments are conducted with [CO₂] fixed at the 2015 level (defined as 2015co2 experiments). Finally, to ensure the availability of picontrol and 2015co2 sensitivity simulations, the emulator calibration is limited to the aforementioned eight GGCMs that participated in ISIMIP3b.

During the preprocessing, all necessary data is sourced from the ISIMIP repository. Consecutive time series (temperatures and precipitations) and static planting and maturity dates from crop calendars are used to derive growing season temperature and accumulated precipitation for each crop. For grid points where the planting date measured by Day of Year (DOY) occurs later than the maturity date, it indicates the crop is planted in one calendar year and harvested in the next. In such cases, the planting year is consistently assigned as the temporal index for all corresponding variables, accounting for the temporal offset in the original GGCM outputs that report consequences of full growing seasons, not of calendar year. GGCMs simulate crop yields across nearly all land areas, even though certain regions do not cultivate all modelled crops, leading to in part unrealistic or invalid data points. Static crop-specific cropland areas (Frieler et al., 2024), irrespective of irrigation conditions, are used as weighting factors to aggregate gridded variables to regional and global scales. This approach excludes grid cells without cropland from the aggregation, thereby improving the representativeness of the results for actual cultivated areas. Detailed information on the crops and corresponding land-use categories is provided in Table S2.

2.3 Updated regional climate

In the OSCAR-crop v1.0, regional growing season climate is derived from the regional climate, which can be provided by OSCAR. To ensure the coherence between the emulator and OSCAR, we update the existing regional climate parameters in OSCAR and further establish relationships between regional and growing season climate variables in the crop emulator. This update is necessitated by two factors: (1) OSCAR v3 utilizes regional climate variables calibrated on CMIP5 data, whereas the Phase 3 GGCMs employ CMIP6 climate data which exhibit higher [CO₂] and warming levels (Jägermeyr et al., 2021; Tebaldi et al., 2021); and (2) the crop emulator operates at a sub-national resolution, creating a scale mismatch with OSCAR v3's broader regional scale (Gasser et al., 2017). Therefore, recalibrating the regional climate is essential for harmonizing the input data across models and resolutions.

A pattern scaling approach (Herger et al., 2015) is employed to derive the relationship between annual global (subscript _G) climate (temperature and precipitation, T_G and P_G) and annual regional climate ($T_{\mathrm{L}}^i$ and $P_{\mathrm{L}}^i$) at a sub-national scale across climate scenarios (Gasser et al., 2017). Herein, the baseline climate is fixed at pre-industrial level, calculated as the average over 1850–2100 from picontrol scenarios. Climate data derived from five ESMs results in five sets of regional climate parameters.

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{L}}^i={\mathit{\omega }}_{T_{\mathrm{L}}}^i\mathrm{\Delta }{T}_{\mathrm{G}}\text{(2)}& {\displaystyle }\mathrm{\Delta }{P}_{\mathrm{L}}^i={\mathit{\omega }}_{P_{\mathrm{L}}}^i\mathrm{\Delta }{P}_{\mathrm{G}}\end{array}$$

Then, the relationship between region-specific (superscript ⁱ) land surface (subscript _L) climate and crop-specific (superscript ^c, including mai, ri1, ri2, soy, swh, and wwh) growing season (subscript _gs) climate ($T_{\mathrm{gs}}^{i,\mathrm{c}}$ and $P_{\mathrm{gs}}^{i,\mathrm{c}}$) is derived as follows. The growing season variables ($\mathrm{\Delta }{T}_{\mathrm{gs}}^{i,\mathrm{c}}$ and $\mathrm{\Delta }{P}_{\mathrm{gs}}^{i,\mathrm{c}}$) are calculated following the preprocessing steps described in Sect. 2.2:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{gs}}^{i,\mathrm{c}}={\mathit{\omega }}_{T_{\mathrm{gs}}}^{i,\mathrm{c}}\mathrm{\Delta }{T}_{\mathrm{L}}^i\text{(4)}& {\displaystyle }\mathrm{\Delta }{P}_{\mathrm{gs}}^{i,\mathrm{c}}={\mathit{\omega }}_{P_{\mathrm{gs}}}^{i,\mathrm{c}}\mathrm{\Delta }{P}_{\mathrm{L}}^i\end{array}$$

These parameters are calibrated based on ISIMIP3b climate scenarios across historical and future periods. The chosen concatenation scheme for the fitting process is to merge all historical and future data directly. We demonstrate the limited impact of this choice in the Supplementary Information (SI). Our investigation considers two schemes for data concatenation. First, a direct merge of all historical and future datasets; second, concatenating each future period specifically with its historical counterpart before merging. As illustrated in Figs. S3–S6, the fitting results show only minor sensitivity, indicating that the fitted parameters are largely unaffected by the choice of scheme. Given this minimal impact, the direct merging scheme is adopted for all calibration procedures due to its simplicity.

3 Model description and calibration

3.1 Response functions

The crop yield responses (R^c) are calculated by dividing crop yields (YD^c) under any given climate and management scenarios by a reference crop yield under a baseline scenario (${\mathrm{YD}}_{\mathrm{0}}^{\mathrm{c}}$), as expressed by Eq. (5):

$$\begin{array}{}\text{(5)}& R^{\mathrm{c}}={\mathrm{YD}}^{\mathrm{c}}/{\mathrm{YD}}_{\mathrm{0}}^{\mathrm{c}}\end{array}$$

The baseline scenario represents pre-industrial temperature and precipitation conditions (as in OSCAR) combined with 2015 levels of [CO₂] and direct human forcing (as defined in the ISIMIP3 experiment design). The total crop yield response can be further broken down into responses to four drivers: [CO₂] ($R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{s}}$), growing season temperature ($R_{T_{\mathrm{gs}}}^{\mathrm{c}}$), accumulated precipitation ($R_{\mathrm{Pgs}}^{\mathrm{c}}$), and nitrogen input ($R_{\mathrm{N}}^{\mathrm{c}}$):

$$\begin{array}{}\text{(6)}& R^{\mathrm{c}}=R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}\times R_{T_{\mathrm{gs}}}^{\mathrm{c}}\times R_{\mathrm{Pgs}}^{\mathrm{c}}\times R_{\mathrm{N}}^{\mathrm{c}}\end{array}$$

It should be noted that all the responses are positive by definition. A driver-specific crop yield response with a value greater than 1 indicates a net positive yield impact, while a value less than 1 indicates a net negative yield impact.

Functional forms for crop yield responses differ across crops, irrigation systems, drivers, and regions, as the fundamental purpose of the emulator is to replicate the outputs of process-based crop models with the highest possible fidelity. To fit the specific response function for each crop-irrigation-driver-region combination, we evaluate seven distinct functional forms, including linear and non-linear regression models, to capture monotonic or unimodal relationships between independent and dependent variables. Under reference conditions, the calculated crop yield responses using these functional forms are intrinsically constrained such that the resulting crop yield response equals 1. Detailed information about functional forms and their mathematical properties is given in Table S3. The functional forms represent four characteristic curve shapes: linear (1), exponential (2), peak-and-decline (3, 4, 6), and level-off (5, 7). Functions (2) through (7) are designed to produce non-negative outputs. When function (1) is employed, any negative yield response values are set to zero.

The calibration of the crop emulator represents a substantial advancement in both functionality and regional specificity. Unlike other modules of the current OSCAR model, which apply uniform functional forms across regions and other module-specific dimensions (optimizing for computational efficiency at the expense of regional relevance), our emulator employs different functional form combinations for each region. This approach enhances the accuracy of regional impact representation, enabling a more nuanced and reliable simulation of crop yield responses under diverse environmental and management conditions.

3.2 Emulator development

3.2.1 Initial quality control

The pre-processed data undergoes initial quality control to remove anomalies or extreme values that could compromise model fitting. After the preprocessing of regional aggregation, the ratios between regional crop yields and their corresponding 5-year rolling averages are calculated before fitting the yield. Data points exhibiting ratios exceeding one order of magnitude (either > 10 or < 0.1) are identified as anomalies in this study and excluded from the training data (Choi et al., 2021). Additionally, all yield values below 1.0 × 10⁻⁶ tDM ha⁻¹ are removed from the dataset to enhance fitting stability. 5-year rolling averages of input and output variables are utilized for subsequent GGCM-based calibration, filtering out short-term data oscillation, as our primary goal is to capture long-term trends and not interannual variability. In the following sections involving GGCMs-related data, original data refers to the regionally aggregated and quality-controlled data, unless stated otherwise.

3.2.2 Parameter limits

To ensure physical validity, parameter limits are integrated into the fitting process of functional forms. Loose upper and lower bounds are imposed on all parameters during the fitting process, informed by prior knowledge regarding the functional form characteristics and the physical interpretation of input variables. For example, in the bell-shaped function (3) (Table S3), the range of parameter β, representing the optimal input level for crop yields, is limited based on established domain knowledge, such as the temperature range for crop growth (Table S4).

3.2.3 Performance metrics

Two metrics, R-squared (R²) and Bayesian Information Criterion (BIC), offering complementary insights, are used to evaluate the emulator's ability to replicate GGCMI Phase 3 simulations during calibration. R², quantifying the variance explained by the emulator, verifies the emulator's ability to reproduce the behavior of complex crop models. Functional forms yielding R² values outside the 0 to 1 range, indicative of inappropriate model fitting, are excluded from the final selection. BIC, a widely used statistical criterion for model selection, balances goodness of fit with model complexity by introducing a penalty term proportional to the number of parameters in the model (Hansen, 2008). This penalty discourages overfitting, ensuring that the selected model not only captures the underlying relationships between variables but also generalizes well to unseen data. By minimizing the BIC scores, we identify functional forms that achieve a trade-off between accuracy and parsimony, ensuring the robustness and computational efficiency of the crop emulator for large-scale applications. The final selection of the optimal functional form for each crop-irrigation-driver-region combination is based on the lowest BIC score identified among all evaluated combinations. The model selection workflow is detailed in Fig. 2.

Figure 2 Flowchart of model development procedure.

The final selected emulator is further validated against the original GGCMs using an additional metric, the relative root mean square error (RRMSE), which quantifies the average magnitude of errors normalized by the mean of the reference values. This normalization enables comparison across crops, regions, and irrigation setups with different yield magnitudes and emphasizes proportional rather than absolute deviations. RRMSE is widely used for model validation in agriculture and climate sciences (Falconnier et al., 2020; Yang et al., 2025; Zheng et al., 2025).

3.3 Main drivers

3.3.1 Atmospheric CO₂ concentrations

[CO₂] is an essential driver of crop yield change through the carbon fertilization effect. Many experiments including Free Air CO₂ Enrichment (FACE) and Open-Top Chamber (OTC) experiments have proved positive crop yield response to elevated [CO₂] (Abebe et al., 2016; Ainsworth and Long, 2021; Bunce, 2016; Cai et al., 2016; Ruiz-Vera et al., 2015). These experiments are conducted under discrete [CO₂] levels (normally one ambient [CO₂] level and one or two elevated [CO₂] levels) using different crop cultivars, mostly in temperate regions with varying soil characteristics, elevation, and regional climate, making it difficult to isolate the impacts of [CO₂] across experiments and establish the relationship between [CO₂] and crop yield response. Additionally, the experimental elevated [CO₂] typically ranges from 550 to 700 ppm, significantly lower than the peak [CO₂] assumed in high-emission scenarios like SSP585 (Toreti et al., 2020). Under the new phase of GGCMI, exploring sensitivities of crop yields to [CO₂] to account for its large uncertainties became one of the main goals (Jägermeyr et al., 2021).

In the crop emulator, ΔCO₂ is the difference between [CO₂] and its reference value at the 2015 level (400 ppm) to ensure consistency with the ISIMIP3b sensitivity experiments, as expressed in Eq. (7). This value is further scaled by reference [CO₂] for model input. Notably, the difference in pre-industrial [CO₂] between ISIMIP3b (285 ppm) and OSCAR (278 ppm) arises from the use of different baseline years – 1850 in ISIMIP3b versus 1750 in OSCAR. Maximum [CO₂] in ISIMIP3b experiments reaches 397 ppm (historical), 474 ppm (SSP126), 867 ppm (SSP370), and 1135 ppm (SSP585).

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}={\mathrm{CO}}_{\mathrm{2}}-{\mathrm{CO}}_{{\mathrm{2}}_{\mathrm{2015}}}\end{array}$$

Then, $R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}$ is quantified as the ratios of crop yields between default sensitivity scenarios and corresponding 2015co2 scenarios per region, so that the impact from temperature and precipitation can be offset according to Eq. (6). Figure S7 presents the cumulative probability distribution of regional CO₂-induced yield responses across all scenarios and regions, disaggregated by GGCM, irrigation and crop type. The values of 95th percentiles of yield responses across crops and GGCMs remain below 3 (typically value of the response < 2). While extreme yield responses (defined as values > 10) account for only 0.02 % of the data, their potential to skew global aggregated results necessitates further investigations. In addition to extreme values inherited from the original data, the emulated yield response can exceed the threshold of 10 even when the original data remains below this level. We therefore consider extreme-value regions arising from both the original data and the best-fit functions. Table S5 provides a detailed breakdown of the specific GGCM-crop-irrigation combinations with regional occurrences associated with these extremes. Only EPIC-IIASA do not produce extreme regional yield responses. Extreme-value occurrences are substantially fewer for firr crops than for noirr crops. Among noirr crops, CYGMA1p74 (swh) and LPJmL (ri1 and ri2) show the highest regional occurrences of extremes. A sensitivity analysis of extreme-value regions on global aggregated result is given in Sect. S1 in the Supplement.

Figure 3 Kernel density estimation of global crop yield responses to [CO₂] change (ppm) across eight GGCMs under historical, SSP126, SSP370, and SSP585, with extreme-value regions masked. Black lines indicate aggregated results from the crop emulator. Green shades and orange shades represent firr and noirr conditions, separately. Darker shades suggest higher density.

We compare original and emulated global yield responses both with (Fig. 3) and without (Fig. S8) masking these extreme-value regions. Kernel density distributions of ΔCO₂ and $R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}$ (shaded areas), derived from original simulations, are shown alongside emulator-based global results (black lines). Reflecting similarities between Figs. 3 and S8, most data cluster around 500 ppm, indicating a strong influence of higher [CO₂] in SSP370 and SSP585 on yield responses. CO₂ fertilization exerts a more pronounced impact on noirr crop yields than firr crop yields, aligning with the water-saving effects of CO₂ fertilization besides higher photosynthesis rates in C3 crops (Deryng et al., 2016; Elliott et al., 2014), though representations of CO₂ fertilization effects vary significantly among GGCMs (Jägermeyr et al., 2021; Rosenzweig et al., 2014). Among the crops herein, maize, a C4 crop, exhibits a relatively small yield response to CO₂. While a general positive trend of CO₂-yield correlation exists, peak-and-decline patterns occur for mai in CYGMA1p74 (both irrigation conditions) and for mai and swh in PEPIC under firr, with SIMPLACE-LINTUL5 simulating a significant decrease in mai yield responses at high [CO₂] under firr. Masking extreme-value regions improves the stability of global yield response emulation at high [CO₂] for swh in CYGMA1p74 and ri1 in LPJmL than unmasked results. Although the emulation's performance at elevated [CO₂] in extreme-value regions is uncertain, we retain the calibrated parameters for these regions in the crop emulator because they reflect the original GGCM simulations. Users of the crop emulator can choose whether to include or disregard results from these regions based on their needs. To mitigate potential bias, this study excludes these regions from regional aggregation during model validation. The strong agreement between kernel density distributions and emulated curves demonstrates the crop emulator's effectiveness in accurately reproducing GGCM behavior.

CO₂ impact on crop yields varies regionally as it is subject to the overall climate regime and crop management (Deryng et al., 2016; McGrath and Lobell, 2013). Figure S9 presents emulated regional crop yield responses to ΔCO₂ across GGCMs under SSP585. The estimation of end-of-century (2069–2099) crop yield responses is performed by computing the mean yield response under future scenarios. End-of-century standard deviations indicate that emulated yield responses for noirr crops demonstrate greater uncertainty than those for firr crops. Of all emulated crops, soy exhibits the largest range of yield responses across GGCMs and regions, indicating larger inter-model uncertainty, while mai shows the narrowest range, suggesting better consistency. Despite the general trend of yield gains due to CO₂ fertilization among all crops ($R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}>\mathrm{1})$, the standard deviation of net yield responses ($R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}-\mathrm{1})$ is greater than their average for certain crop-region combinations. For instance, mai yields under firr in all the selected regions, and soy yields under noirr in Brazil and India display net yield response ranges with mixed signs, indicative of potential yield reductions under high [CO₂].

3.3.2 Climate

Climate variability accounts for a significant percentage of annual crop yield variance (Anderson et al., 2019; Ray et al., 2015), and the ongoing warming trend exacerbates the yield fluctuations (Ostberg et al., 2018). Field warming experiments generally demonstrate negative yield responses for maize, rice, soybean and wheat under elevated temperatures (Wang et al., 2020; Zhao et al., 2016). Furthermore, the couplings between rising temperature and decreased moisture intensify crop yield sensitivities (Lesk et al., 2021). While extreme weather events pose significant threats to crop production, their impacts and associated crop yield responses are underestimated by climate and crop models (Kornhuber et al., 2023), a factor not explored in this study. In water-limited regions, irrigation is an intensification strategy to enhance crop production under various levels of warming. However, its potential in buffering climate change is constrained by water availability, thereby limiting its agricultural benefits (Elliott et al., 2014; Li and Troy, 2018; Minoli et al., 2019; Wang et al., 2021).

The crop emulator utilizes temperature differences ($\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}$) and relative precipitation changes ($\mathrm{\Delta }{P}_{\mathrm{gs}}^{\mathrm{c}}/P_{\mathrm{gs},\mathrm{pi}}^{\mathrm{c}}$) as input variables. Reference growing season temperature ($T_{\mathrm{gs},\mathrm{pi}}^{\mathrm{c}}$) and precipitation ($P_{\mathrm{gs},\mathrm{pi}}^{\mathrm{c}}$) are obtained from picontrol climate scenarios.

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}=T_{\mathrm{gs}}^{\mathrm{c}}-T_{\mathrm{gs},\mathrm{pi}}^{\mathrm{c}}\text{(9)}& {\displaystyle }\mathrm{\Delta }{P}_{\mathrm{gs}}^{\mathrm{c}}=P_{\mathrm{gs}}^{\mathrm{c}}-P_{\mathrm{gs},\mathrm{pi}}^{\mathrm{c}}\end{array}$$

For firr crops, where water availability is non-limiting, the water stress term $R_{\mathrm{Pgs}}^{\mathrm{c}}$ in Eq. (6) is excluded (taken equal to 1). Given that nitrogen input is fixed at the 2015 level, meaning $R_{\mathrm{N}}^{\mathrm{c}}$ also equals to 1 in current ISIMIP3 simulations, $R_{\mathrm{Tgs}}^{\mathrm{c}}$ is equivalent to the ratio of composite R^c to the fitted $R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}$ for firr crops with $R_{\mathrm{Pgs}}^{\mathrm{c}}=\mathrm{1}$. For noirr crops, where the water stress term is not constant, the functions of growing season temperature and precipitation are fitted simultaneously through multivariate regression, as experiments isolating the individual impacts of temperature and precipitation are lacking. The combined effect of $R_{\mathrm{Tgs}}^{\mathrm{c}}$ and $R_{\mathrm{Pgs}}^{\mathrm{c}}$ is calculated as $R^{\mathrm{c}}/R_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{c}}$. Seven functional forms (Table S3) are applied to each driver, resulting in 49 and 343 possible functional form combinations for firr and noirr crops, respectively, per region. The optimal combinations of functional forms are subsequently determined based on the selection criteria described above.

Figure 4 Global crop yield responses to growing season temperature change (K) across eight GGCMs (extreme-value regions masked), showing firr (green solid lines) and noirr (orange dashed lines) crop yield response sensitivity to growing season temperature change under SSP585, and RRMSE (%) between firr and noirr crops across GGCM-crop combinations. Sensitivity values for firr and noirr crops are shown in green and orange, respectively, with their RRMSE values reported in the upper-right corner.

Figure 4 displays the emulated global crop yield responses to $\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}$ under SSP585 which shows the widest temperature change among scenarios. Across all GGCMs, rising $\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}$, caused by increasing land temperature, generally drives yield declines in the emulation. However, crop yields exhibit peak-and-decline trends as growing season temperature increases for some GGCM-crop combinations such as all SIMPLACE-LINTUL5-based crops and PROMET-based wwh. Furthermore, the sensitivity of yield response is quantified as the ratio of yield response differences between 2100 and 2015 to corresponding growing season temperature differences. Under temperature change in SSP585, CYGMA1p74-based emulations show the strongest yield response sensitivities for all crops, except for swh. Among all crops, soy is the most sensitive to rising temperature, while wwh demonstrates more consistent and less sensitive responses. Despite a moderate average response, mai shows more divergent response sensitivities across GGCMs. Globally, firr and noirr crops exhibit generally comparable yield responses induced by $\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}$ across most GGCM-crop combinations, especially for mai in LPJmL and swh in PROMET, as indicated by the smallest RRMSE. More significant divergences emerge for ri1 in PROMET, soy in SIMPLACE-LINTUL5, and swh in CYGMA1p74 and PEPIC.

Regional emulations of crop yield responses to $\mathrm{\Delta }{T}_{\mathrm{gs}}^{\mathrm{c}}$ are shown in Fig. S10. End-of-century crop yield responses across the majority of GGCM-crop-irrigation-region combinations suggest a temperature-induced decrease in mai, ri1, ri2, soy, and swh yields across major producing regions, compared to the pre-industrial reference climate. In contrast, potential yield gains are projected for wwh in Russia, predominantly influenced by high emulation results derived from PROMET. Among all crops, end-of-century soy yield is the most negatively impacted by rising temperature, especially in Brazil, the USA, and India, while wwh yield is more resistant to temperature change in most major producing regions, except for India.

Figure 5 Global crop yield responses to relative precipitation change (%) across eight GGCMs under SSP585 (extreme-value regions masked).

The aggregated emulations of global crop yield responses to relative precipitation changes under SSP585 show considerable variability, as depicted in Fig. 5. Most GGCM-crop combinations indicate positive yield responses to increased precipitation. Notably, there is a significant increasing trend of swh yield responses across all GGCM-based emulations. However, global soy yield responses emulated from CYGMA1p74 decrease with enhanced precipitation, consistent with the findings that excessive water reduces legume yields (Iizumi et al., 2024). Additionally, global mai and soy yield responses show substantial inter- and intra-GGCM uncertainties. These patterns underscore the necessity of representing the spatially heterogeneous nature of precipitation changes and corresponding impacts.

Figure S11 illustrates the emulated regional noirr yield responses to relative growing season precipitation change. Under SSP585, global end-of-century swh yield exhibits the greatest response to precipitation change, followed by mai. At the regional level, mai and swh in EU27 show significant negative impacts from precipitation, whereas mai and soy in Argentina, and swh and wwh in China benefit from precipitation increase. Notably, SIMPLACE-LINTUL5-based emulations consistently predict the most pronounced yield responses to precipitation for mai, swh and wwh across most major production regions.

3.3.3 Nitrogen input

Nitrogen fertilizer significantly influences crop yields, particularly in nitrogen-limited regions like Africa (Falconnier et al., 2020; Olin et al., 2015). Nitrogen management, balancing use efficiency and environmental impact, has become a central focus in the development of sustainable agricultural practices (Springmann et al., 2018; Wang et al., 2019). While the interaction between the nitrogen cycle and climate affects the nitrogen uptake by crops (Liu et al., 2024; Wang et al., 2025), GGCM-based simulations exploring crop yield responses to varying nitrogen fertilizer applications remain insufficient to support robust emulation. During GGCM Phase 1, historical crop yields were simulated under three harmonization levels with model-specific assumptions about nitrogen stress (Elliott et al., 2015; Müller et al., 2019). GGCMI Phase 2 simulations evaluated crop yield responses to three discrete nitrogen fertilizer input levels; however, the limited sample size precludes reliable model fitting (Franke et al., 2020a, b). In GGCM Phase 3, the impact of nitrogen fertilizer is not considered in the available simulations, as fixed nitrogen input levels according to present spatial patterns are applied across all historical and future scenarios.

To facilitate nitrogen impact assessment in the crop emulator, data from long-term field experiments are employed to quantify crop yield responses to nitrogen inputs following van Grinsven et al. (2022), where the relationship between nitrogen inputs and crop yield response is described by a quadratic function for cereal crops, assuming constant soil nitrogen supply during the experimental periods. As for nitrogen-fixing soybeans, biological nitrogen fixation (BNF) exhibits a negative correlation with fertilization rate (Ma et al., 2022; Salvagiotti et al., 2008), leading to large uncertainty in the prediction of soybean yield response to nitrogen (Luiz et al., 2023; Pannecoucque et al., 2022). Furthermore, the variability in soybean yields attributed to nitrogen management is marginal compared to climatic and irrigation factors (Mourtzinis et al., 2018). Consistent with this, an existing GGCMI Phase 2 crop emulator indicates a very limited GGCM-based nitrogen response for soybean (Franke et al., 2020a; Liu et al., 2023). Therefore, nitrogen input is not considered a yield determinant for soybean in our crop emulator for now, resulting in a constant yield response ($R_{\mathrm{N}}^{\mathrm{soy}}=\mathrm{1}$).

For cereal crops in the emulator, total nitrogen input comprises fertilizer application, BNF, and atmospheric deposition (Zhang et al., 2015). Nitrogen fertilization rates are spatially and temporally variable throughout the historical period. BNF for non-symbiotic cereal crops is assumed to be crop- and region-specific and constant over time (Ladha et al., 2022). Also, nitrogen deposition is treated as static, despite observed temporal variations (Ackerman et al., 2019; Yu et al., 2019).

$$\begin{array}{}\text{(10)}& N_{\mathrm{av}}=N_{\mathrm{fer}}+N_{\mathrm{man}}+N_{\mathrm{bnf}}+N_{\mathrm{dep}}\end{array}$$

Crop-specific response functions are refitted using Michaelis-Menten, Mitscherlich, and George models (Table S6), replacing the generic second-order polynomial form utilized in the original study to ensure non-negative yield responses. The optimal functional form is subsequently selected based on the R² and BIC criteria discussed previously. Due to the limited size of sampling data, the response function for the nitrogen driver is not region-specific. To ensure consistent estimation of total nitrogen input, the extrapolated left x-intercept of the second-order polynomial function is used to approximate soil nitrogen supply during the fitting of nitrogen yield responses using the three alternative models. The nitrogen inputs are normalized by a scaling factor of 100 kgN ha⁻¹ ($N_{\mathrm{av}}/\mathrm{100}$) to ensure physical interpretability of the exponent in the George model. As illustrated in Fig. 6, both 2nd order and George models exhibit peak-and-decline yield response. Michaelis-Menten and Mitscherlich models level off after nitrogen input exceeding 200 kgN ha⁻¹, aligning with the saturation effect of nitrogen input for GGCMI Phase 2 emulators (Franke et al., 2020a). Similarly, simulated global wheat yield responses to increasing nitrogen fertilizer application rates also saturate at approximately 200 kgN ha⁻¹ for wheat (Martre et al., 2024).

Figure 6 Crop-specific yield responses to nitrogen inputs (100 kgN ha⁻¹). Maize and wheat data are from van Grinsven et al. (2022) Table S2; rice data compiled from sources in Note 8.

While there might be negative impacts of nitrogen fertilizer on crop yields (van Grinsven et al., 2022; Maaz et al., 2021), current GGCMs don't incorporate functions resulting in yield penalties at high nitrogen inputs. Moreover, excessive nitrogen application rate above 450 kgN ha⁻¹ results in zero yield for wheat in the 2nd order and George models, which is not typically observed in experiments or real-life applications (Yokamo et al., 2023). Therefore, Mitscherlich function with better performance is selected to represent the nitrogen-yield response in the crop emulator.

4 Emulator validation

4.1 In-sample validation

To examine the validity of the crop emulator, the RRMSEs between our emulations and the original ISIMIP3b GGCM simulations across all climate scenarios are calculated at different regional levels. As depicted in Fig. 7, RRMSEs for global crop yields across all GGCM-crop-irrigation combinations remain below 12 %, with most values clustered below 5 %. Under firr, mai exhibits the widest range of RRMSE across GGCMs, whereas soy and swh show tighter distributions. Under noirr, swh displays the largest spread and the highest RRMSE derived from SIMPLACE-LINTUL-based results, while the RRMSE distribution for ri1 is the most compact. The lowest errors are consistently associated with LDNDC-based mai, ri1, ri2, soy and swh yields under firr. Overall, our crop emulator reproduces the GGCM outputs with high fidelity on a global scale.

Figure 7 RRMSE (%) between global area-weighted emulated and original crop yields across crop-GGCM combinations under historical, SSP126, SSP370, and SSP585 scenarios. For each crop, the markers along the horizontal axis follow the alphabetic order of the eight GGCMs.

Moreover, the RRMSE for the six producing regions is visualized as heatmaps in Fig. 8. Most GGCM-irrigation-crop-region combinations exhibit low RRMSE values (< 10 %) between the emulated and original simulation results, especially for ri1 yields which displays good consistency across GGCMs, irrigation regimes, and regions. Regionally, the RRMSE values for the rest of the world (R6) are consistently low. Under noirr, swh and wwh yields in India, as well as ri2 yields in Bangladesh demonstrate relatively large deviations from original GGCMs. Under firr, wwh yields in India exhibit moderate RRMSEs. Notably, CYGMA1p74-based swh yields in India under noirr shows the highest RRMSE, reaching 41 %. In addition, SIMPLACE-LINTUL5-based emulations for mai, swh, and wwh under noirr show comparably high RRMSE values. Overall, our crop emulator demonstrates robust performance over the most important producing regions.

Figure 8 RRMSE (%) between emulated and original crop yields in six producing regions (R1 to R6 from Table S1) for each crop-GGCM combination under historical, SSP126, SSP370, and SSP585 scenarios.

To provide a more detailed assessment of the crop emulator's performance at its native spatial resolution, Figs. S12–S19 present sub-national, end-of-century crop yield differences between our emulations and original simulations under SSP585. These figures are accompanied by the global weighted absolute and relative yield differences, the latter given in percentages. Under SSP585, the extreme changes in [CO₂], temperature, and precipitation lead to correspondingly large variations in crop yield responses. At the crop emulator's finest regional scale, absolute crop yield differences are less than 0.5 tDM ha⁻¹ for the majority of GGCM-crop-irrigation combinations, demonstrating the emulator's capability to reproduce the results of GGCMs at a sub-national level. When aggregated globally, the largest absolute global yield difference occurs in SIMPLACE-LINTUL5 for swh yield under noirr, resulting in a relative yield difference of less than 8 % between our emulator and the original GGCM (Fig. S19), while ISAM for mai yield under noirr shows the smallest relative yield difference (Fig. S14). The low percentage of yield differences relative to original yields suggests that the discrepancies between emulated and original end-of-century yields are, for the most part, negligible. Overall, the emulation results at both sub-national and aggregated levels show good consistency with the original GGCMs' simulations in terms of general trends and absolute values.

These results also illustrate the interaction effects inherent in our multiplicative model structure. The SSP585 scenario represents a high-stress regime where elevated CO₂, high temperatures, and potential drought conditions co-occur – precisely the conditions where interactions between drivers are expected to be largest. The good agreement between emulator and GGCM outputs under this extreme scenario indicates that the implicit interactions arising from our nonlinear component functions capture the coupled responses present in the original process-based models.

4.2 Out-of-sample validation

To further validate our crop emulator, historical simulations are performed using observation-driven ISIMIP3a GSWP3-W5E5 climate data as inputs. Daily climate variables are preprocessed into annual resolution and then offset by subtracting their averages over the 1901–1930 period (approximating pre-industrial baseline climate). The pattern scaling method (Sect. 2.2) is then applied to derive regional and growing season climate variables within the crop emulator. The resulting outputs, averaged over the period of 1980–2010, are compared against original ISIMIP3a historical simulations. Because the emulator is calibrated using simulations from multiple GGCMs simultaneously, it is designed to reproduce the collective response space of the GGCM ensemble rather than the behavior of any individual model. The out-of-sample validation therefore evaluates emulator performance under independent climate forcings while remaining within the range of crop responses represented by the GGCM ensemble.

Global maps in Fig. 9 show the ensemble means of the yield differences between our emulations and GGCM simulations, which are generally within ±0.5 tDM ha⁻¹ among most national or sub-national regions. Among top five producing regions, the most significant underestimation of multi-model mean occurs for mai yields under noirr in Argentina with an absolute yield difference of −0.6 tDM ha⁻¹. The violin plots in Fig. 9 reveal that global SIMPLACE-LINTUL5-based crop yields exhibit the largest relative deviations from the original GGCM under both firr and noirr, except for ri1 and ri1 which SIMPLACE-LINTUL5 does not simulate. For ri1 and ri2, the PROMET-based and CYGMA1p75-based results show the most significant deviations from the original GGCMs under firr and noirr, respectively. Notably, the crop emulator demonstrates the highest consistency with GGCMs for global wwh yields under noirr, indicated by the smallest range of global relative yield difference in the violin plot, while the largest discrepancy is observed for global ri2 yields under firr. Generally, the multi-model means show small differences between the crop emulator and original GGCMs during the recent decades.

Figure 9 Yield differences between crop emulator and original ISIMIP3a historical simulations shown as: global maps of GGCM-averaged sub-national absolute differences (tDM ha⁻¹) and violin plots of relative differences (%) for global average values across crop-GGCM combinations. For each crop in the violin plots, the markers along the horizontal axis follow the alphabetic order of the eight GGCMs.

Individual GGCM-emulator yield discrepancies are displayed in Figs. S20–S27. Our emulator exhibits the largest absolute underestimation for SIMPLACE-LINTUL5 mai yields under firr, where the percentage of global yield difference is approximately −13 % (Fig. S27). Nevertheless, the most significant relative yield difference is observed for PROMET ri2 yields under firr, where the relative percentage difference is around 14 % (Fig. S26). Compared with original GGCMs, the crop emulator tends to overestimate CYGMA1p74 crop yields in the recent historical period, with a less than 10 % increase in the global crop yield estimations (Fig. S20). On the contrary, SIMPLACE-LINTUL5 regional emulations are generally lower than the original results (Fig. S27). The crop emulator demonstrates high performance for EPIC-IIASA, LDNDC, and PEPIC emulations, with absolute yield differences being less than 5 % of the GGCMs' results. A comparison of the multi-model ensemble average with individual GGCM outputs demonstrates both the robustness of ensemble-based projections and the emulator's fidelity in replicating complex model behavior.

5 First results

5.1 Comparison with field experiments

To ground the emulator's projections in real-world experimental evidence, distributions of crop yield sensitivities derived from FACE, OTC, and field warming experiments are compared with diagnostic outputs from the crop emulator. This sensitivity analysis is conducted at global scale, with regional outputs from the emulator aggregated before comparison. For CO₂ fertilization effects, crop yield sensitivities to elevated [CO₂] are quantified by first computing the relative yield difference between a historical baseline period (1991–2010) and a future period under the SSP585 high-emissions scenario, characterized by the experiments' elevated [CO₂] levels. This relative difference is then normalized by the [CO₂] increment, thereby aligning the emulator-based sensitivities with the experimental design of FACE and OTC studies in which ambient [CO₂] typically ranges from 320–420 ppm and elevated levels span 500–800 ppm. A 20-year analysis window under SSP585 is dynamically selected to match the [CO₂] enrichment levels observed in experiments. For example, if experimental studies for maize report an average enrichment of +200 ppm, the emulator's future period is chosen as the 20-year interval in which SSP585 [CO₂] levels are closest to the historical baseline plus 200 ppm. A similar approach is applied to temperature sensitivity. Experimental yield sensitivities to warming are compared against emulator-based sensitivities, calculated as the relative yield change per degree Kelvin increase in growing season temperature.

To generate the probabilistic distribution of yield sensitivities, the emulator is run with a Monte Carlo (MC) ensemble of 1000 parameter configurations. Input datasets are consistent with those of ISIMIP3b. Parameters are randomly drawn from five regional and growing-season parameter sets and eight crop yield response parameter sets. To extend parametric uncertainty, each sampled parameter value is perturbed by an additional noise covering a 10 % relative standard deviation, thereby capturing a wider range of plausible yield responses.

Figure 10 Frequency (patches) and corresponding fitted normal distributions (lines) of global crop yield sensitivities to elevated [CO₂] and growing season temperature from experiments (dotted blue lines) and simulations (solid lines) from the crop emulator.

Figure 10 presents the distributions of global crop yield sensitivities to elevated [CO₂] and growing season temperature. Overall, the emulator's yield sensitivity estimates demonstrate good agreement with experimental results in terms of the average values. Compared to field experiments, emulator-derived sensitivities exhibit narrower distributions, suggesting the GGCMs may not cover the full range of possibilities measured in the field. Field experiments reveal greater uncertainty in maize crop yield sensitivities to [CO₂], with the emulator-based sensitivities aligning with the majority of experimental data. While the experimental data employed does not distinguish between the two rice growing seasons or the spring and winter wheat types, higher sensitivity to [CO₂] for ri2 than ri1 is identified based on the crop emulator, whereas the sensitivity to growing season temperatures is comparable for both growing seasons. For the two wheat types, swh displays a higher sensitivity to [CO₂] than wwh, and its negative response to growing season temperature for swh is also more pronounced than that of wwh. The substantial consistency observed between the emulated and experimentally derived crop yield sensitivities provides compelling evidence for the robust performance of the emulator in representing real-world crop yield behavior under changing climatic conditions, although this could be expected as GGCMs are calibrated on this or similar observational data.

5.2 Comparison with FAOSTAT

Historical simulations incorporating dynamic climate forcing and both time-varying or fixed 2015-level nitrogen inputs are conducted to enhance our crop emulator's real-life applicability. Temperature and precipitation inputs are taken from 20CRv3-ERA5 (Compo et al., 2011; Hersbach et al., 2020; Slivinski et al., 2019) and CRU-TS4.09 (Harris et al., 2020) for the period 1901–2021, with the 1901–1930 average climate used as the pre-industrial baseline to calculate temperature and precipitation changes. Crop-type-specific nitrogen fertilizer application rates from LUH2 (Hurtt et al., 2020) and crop-specific nitrogen fertilizer application rates from GGCMI Phase 3 (Heinke et al., 2021) are employed to better represent regional and crop-type disparities in nitrogen inputs. Since both datasets lack data beyond 2015, nitrogen inputs for 2015–2021 are estimated using 5-year rolling totals from FAO, scaled by the nitrogen application ratios between FAO and the two datasets over the 2010–2014 period. For direct comparison with FAO statistics, simulated yields for ri1 and ri2 are aggregated based on their respective harvested areas (Jägermeyr et al., 2021; Laborte et al., 2017; Waha et al., 2020), with a similar aggregation approach applied to swh and wwh (Becker-Reshef et al., 2023). Furthermore, firr and noirr crop yields are combined based on irrigated fractions of cropland (Frieler et al., 2024). Following this post-processing, the yields of the four major crops are validated against FAO records, using only countries with valid data from both FAO and our emulations to avoid mismatches.

Including nitrogen substantially improves the crop emulator's capability to replicate historical crop yield trends as shown in Fig. 11, where the lines exhibit multi-ensemble mean crop yields, and the shades show the ±1 standard deviation uncertainty range induced by 1000 MC configurations. Without nitrogen inputs, simulated yields exhibit notable discrepancies from FAO data, particularly in capturing temporal trends, highlighting nitrogen's key role in historical yield changes. Maize, rice, and wheat yield trends are well-captured with nitrogen inputs. In some crop-region combinations, such as rice in Indonesia, the emulated crop yields with nitrogen impact are highly consistent with FAO data, indicating a significant improvement in crop yield dynamics after nitrogen impact. Nevertheless, rice yields are well estimated across most major producing regions, while soybean yields are overestimated under the assumption of insignificant nitrogen soybean yield responses as described previously. Stagnant or declining yields in several regions are estimated by the crop emulator since 1990s such as maize in the USA and EU27, in consistent with the findings of yield stagnation in some major crop growing regions (Gerber et al., 2024; Ray et al., 2012), while the FAO data suggest continued increase in yield in the meantime.

Figure 11 Regional crop yields (tDM ha⁻¹) of maize, rice, soybean and wheat simulated by the crop emulator. Blue lines represent simulations with nitrogen impacts, and red lines represent simulations without nitrogen impacts. FAO data are shown as brown dotted lines. The shaded areas indicate ±1 standard deviation of uncertainty arising for multi-model differences embedded in the crop emulator.

Moreover, discrepancies between simulated and FAO yields arise from several factors. While the emulator accounts for nitrogen inputs, it currently does not incorporate other influencing factors such as labor availability, cultivar improvements, adaptation-related crop calendar and growing region shifts, or pesticide/herbicide use, all of which impact yield outcomes (Abramoff et al., 2023; Yang et al., 2024). Also, triple-cropping of rice in South Asia, where rice yields and productions are high, is not represented in current GGCMs thus not included in the crop emulator. Furthermore, the crop emulator's current formulation does not consider the impact of extreme weather events and therefore only represent multi-year trends. This limitation comes from both the modeling paradigm of OSCAR and the limited capacity of GGCMs to accurately simulate these phenomena (Liu et al., 2023). Additional sources of uncertainties encompassing the input datasets, such as fertilizer inputs and irrigated land fractions, also contribute to the discrepancies. Despite these limitations, the crop emulator effectively demonstrates its ability to reproduce some key temporal yield patterns that show good consistency with observations during the historical period.

6 Discussion and conclusion

This study presents the development and evaluation of a sub-national scale crop yield emulator integrated within the compact Earth system model OSCAR. By leveraging outputs from eight state-of-the-art GGCMs driven by five CMIP6 ESMs, the emulator simulates yields of four major crops under varying climatic and management conditions. Through a multi-step calibration process and the statistical selection of empirical functional forms, the emulator captures key biophysical yield responses to [CO₂], growing season temperature, water availability, and nitrogen inputs. The emulator's region-specific functional forms and sub-national resolution offer significant advantages in terms of regional representation, overall performance, and policy relevance.

Performance evaluations against GGCMI Phase 3b (in-sample) and GGCMI Phase 3a (out-of-sample) simulations confirm the emulator's ability to precisely replicate complex crop model behavior. Correlation coefficients across most regions and GGCMs are good, and the emulator effectively reproduces yield distributions at global and regional scales. Importantly, the comparison with field experiments corroborates the emulator's robustness. Comparisons with FAO data further validate the importance of nitrogen representation in the emulator to reproduce observed historical yield patterns, although there remains room to introduce additional driving factors to better reproduce these historical yields.

One of the key strengths of the crop emulator lies in its balance between biophysical accuracy and computational efficiency, making it especially suitable for large-scale scenario analysis, probabilistic impact assessments, and integration into SCMs or integrated assessment models (IAMs). For example, reanalysis of existing scenario databases, such as the Sixth Assessment Report (AR6) scenario database (Kikstra et al., 2022), can provide insights into disparities between projections made by IAMs and GGCMs (Calvin and Fisher-Vanden, 2017). The emulator also enables assessment of negative emissions derived from crop residues utilized for bioenergy under various land-use scenarios (Xu et al., 2022). Coupling or linking the emulator with IAMs can further improve their representation of climate impacts on agriculture (Ruane et al., 2017). The probabilistic framework, based on Monte Carlo sampling, enables ensemble projections that capture uncertainty across multi-model specifications (Gasser et al., 2017; Quilcaille et al., 2023), thereby informing risk assessments and scenario planning (Kephe et al., 2021; Lesk et al., 2022). Moreover, the emulator's regional flexibility and resolution allow for detailed assessments of spatial yield variability and inter-country differences, which are often critical for regional policy and adaptation strategies.

Despite its strengths, the current version of the emulator lacks some socioeconomic drivers such as seed technology (e.g., Hasegawa et al., 2019), pest and disease management (e.g., Sheng et al., 2025), and labor availability (Orlov et al., 2024; Sheng et al., 2025), all of which may influence crop yield dynamics. These factors are closely linked to Agricultural Research and Development (R&D), which is recognized as a key driver of yield growth (Hultgren et al., 2025; Iizumi et al., 2025). Additionally, while the emulator simulates nitrogen impacts using long-term field experiments, uncertainty remains for crops like soybeans that rely on biological nitrogen fixation (Mourtzinis et al., 2018).

While our emulator aims to represent real-world conditions, the fully irrigated (firr) assumption implies irrigation water is supplied whenever required to prevent soil moisture stress. This simplification does not reflect real-world irrigated croplands, where water resources and infrastructure constrain irrigation (Elliott et al., 2014). Moreover, it should be noted that the emulator is suited for applications focusing on long-term climate responses and scenario exploration rather than assessments of extreme events such as droughts and heatwaves (Lesk et al., 2016; Zampieri et al., 2017).

Future work should aim to enhance the crop emulator's versatility and accuracy by (1) incorporating additional agronomic adaptation practices, such as crop switching (Kling et al., 2024); (2) extending the crop emulator to include bioenergy crops to improve its relevance for the IAM community (Ewert et al., 2015); (3) integrating higher-quality inputs including land-use and crop calendar data to improve the parameterization and interpretation (Nóia-Júnior et al., 2025); and (4) implementing crop-specific thresholds for extreme climate events to facilitate more detailed assessments of risk management (Hasegawa et al., 2021; Nóia-Júnior et al., 2025).

In conclusion, the crop emulator developed in this study demonstrates strong capability in reproducing GGCM simulations, offering a computationally efficient yet accurate tool for projecting crop productivity under climate change. Its integration into OSCAR paves the way for comprehensive assessments of future food security and land-based mitigation strategies in Earth system modeling frameworks. Continued refinement and expansion of this emulator will enhance its utility for informing policy and research at global, regional, and sub-national levels.