The statistical emulators of GGCMI phase 2: responses of year-to-year variation of crop yield to CO₂, temperature, water, and nitrogen perturbations

Understanding the impact of climate change on year-to-year variation of crop yield is critical to global food stability and security. While crop model emulators are believed to be lightweight tools to replace the models, few emulators have been developed to capture such interannual variation of crop yield in response to climate variability. In this study, we developed a statistical emulator with a machine learning algorithm to reproduce the response of year-to-year variation of four crop yields to CO₂ (C), temperature (T), water (W), and nitrogen (N) perturbations defined in the Global Gridded Crop Model Intercomparison Project (GGCMI) phase 2. The emulators were able to explain more than 52 % of the variance of simulated yield and performed well in capturing the year-to-year variation of global average and gridded crop yield over current croplands in the baseline. With the changes in CO₂–temperature–water–nitrogen (CTWN) perturbations, the emulators could reproduce the year-to-year variation of crop yield well over most current cropland. The variation of R and the mean absolute error was small under the single CTWN perturbations and dual-factor perturbations. These emulators thus provide statistical response surfaces of yield, including both its mean and interannual variability, to climate factors. They could facilitate spatiotemporal downscaling of crop model simulation, projecting the changes in crop yield variability in the future and serving as a lightweight tool for multi-model ensemble simulation. The emulators enhanced the flexibility of crop yield estimates and expanded the application of large-ensemble simulations of crop yield under climate change.

1 Introduction

The impact of climate change on crop yield is an increasing concern related to global food security (Kinnunen et al., 2020). Two distinct approaches have been used to evaluate the impact of climate change on crop yield, process-based crop models and statistical models. Process-based crop models are reliable tools to project crop yields under future climate change but are computationally expensive (Jones et al., 2017). In contrast, statistical models are lightweight tools that can fit yield response to historical climate change (Li et al., 2019b), but the relationship between climate factors and crop yield is based on the historical climate conditions and their effects on crop yields, which can hardly be used for future projection with new, unprecedented climate conditions. Therefore, it is promising to develop tools that can reduce the expense of computation and increase capacity for flexible future projections (Franke et al., 2020a).

Earlier studies have developed statistical emulators of process-based crop model results to balance the advantages and disadvantages of process-based crop models and statistical models. Those statistical emulators were initially developed with “entire scenario” (simultaneous changes in climate factors) simulation during historical or future periods. Emulators have been developed for process-based crop models, like APSIM (Shahhosseini et al., 2019), GEPIC (Folberth et al., 2019), GWG (Xu et al., 2021), GAZE (Raimondo et al., 2021), and WOFOST (Tartarini et al., 2021), and used to estimate historical crop yield. As an emulator trained by historical simulations cannot project the crop yield in the future, multiple crop model ensemble simulations in future climate scenarios have been used to calibrate emulators (Blanc, 2017, 2020; Blanc and Sultan, 2015; Mistry et al., 2017; Ostberg et al., 2018). However, scenario-based future crop yield projection is not a systematic perturbation of climate factors (Franke et al., 2020a). For instance, scenario-based yield projection can only provide the simulated crop yield driven by simultaneous changes in climate factors. The dependency of temperature and precipitation will be kept in scenarios such that the impact of temperature and precipitation cannot be clearly separated.

An alternative emulation based on a training dataset of “perturbed factors” was introduced, which offers advantages compared to separate effects of crop yield drivers. Perturbed factor emulation was first conducted with site-based crop model simulations, which could estimate the yield across a broad range of CO₂, temperature, and water (Fronzek et al., 2018; Makowski et al., 2015; Pirttioja et al., 2015), but these emulators were limited to the site level. To break the constraints of site-based simulation, the Global Gridded Crop Model Intercomparison (GGCMI) phase 2 provided a simulation dataset across structured CO₂–temperature–water–nitrogen (CTWN) perturbation cubes. This dataset offered two major advantages: it allows for separating the effects of different climatic factors and of nitrogen levels on crop yields, and it allows distinguishing the climatological mean and year-to-year variation of yields (Franke et al., 2020b). The phase 2 dataset was published to support the derivation of crop yield climate change “response surfaces”. Based on the CTWN cubes, a statistical emulator has been developed, providing near-global-coverage multi-model emulators of climatological mean yield projections from the GGCMI phase 2 ensemble by using a regression model with a third-order polynomial basis function (Franke et al., 2020a). Due to the focus on climatological mean yield, the aspect of year-to-year variation of yield under CTWN perturbations has not been fully analyzed or exploited in emulator design.

For climate change risk assessment, interannual yield variability (or the year-to-year variation of yield) is an important metric of yield risk (Liu et al., 2021b) and food supply stability (Liu et al., 2021a) but has been insufficiently addressed in previous studies (Campbell et al., 2016). Large year-to-year variation of crop yield can influence livelihoods of producers, food prices (Hasegawa et al., 2021), and hunger (Janssens et al., 2020) and even lead to political instabilities (Sternberg, 2011). Recently, year-to-year variation has been introduced as a metric for climate change risk to global crop production (Jägermeyr et al., 2021). Developing statistical emulators that can reproduce the year-to-year variation of yield from the CTWN cubes could therefore provide a powerful tool for studies focusing on the risk of climate change impact on yield. In this study, we exclusively aimed to develop statistical emulators to reproduce year-to-year yield variation with GGCMI phase 2 data.

Table 1 GGCMs included in emulation. Each model offers the same set of CTWN simulations across four crops.

^∗ LPJ-GUESS omits maize and rice, and ORCHIDEE-crop omits spring wheat (denoted by ×).

2 Data and methods

2.1 Data

The input and output data for the simulation of global gridded crop yield were obtained from the GGCMI phase 2 dataset, which includes gridded crop yield projections at 0.5^∘ longitudinal–latitudinal resolution for maize, spring wheat, winter wheat, and rice (Franke et al., 2020b). The input data for the process-based simulations in GGCMI phase 2 included data for climate, soil, atmospheric CO₂ concentration, and nitrogen fertilizer application rates. Baseline (1980–2010) climate inputs were used from the AgMIP Modern-Era Retrospective Analysis for Research and Applications (AgMERRA) forcing dataset, including daily maximum and minimum temperatures, precipitation, and solar radiation (Ruane et al., 2015). Systematic perturbations were conducted in each grid cell with seven temperature levels (from −1 to +6 K in 1 K intervals, with +5 K skipped), nine precipitation levels (from −50 % to +30 % in 10 % intervals, with −40 % skipped; the Winf precipitation level presents simulations under fully irrigated conditions), four CO₂ concentration levels (360, 510, 660, and 810 ppm), and three nitrogen levels (10, 60, and 200 kg ha⁻¹). Simulations were repeated for two adaptation strategies, i.e., no adaptation in cultivars (A0) and adaptation by maintaining growing season length (A1). A total of 12 global gridded crop models (GGCMs) were then forced with each of these perturbations of the original reanalysis weather data. We selected 10 of 12 crop models in the GGCMI phase 2 experiment for constructing the emulators. These were APSIM-UGOE, CARAIB, EPIC-IIASA, EPIC-TAMU, GEPIC, LPJ-GUESS, LPJmL, ORCHIDEE-crop, pDSSAT, and PEPIC (Table 1). PROMET and JULES were not included as they used different climate inputs.

The GGCMs used a national and subnational crop calendar based on Sacks et al. (2010), Portmann et al. (2010), and environment-based extrapolations (Elliott et al., 2015). The crop calendar was used to determine the window to calculate the climatic predictors and grid-specific growing season length. The current global harvested area for identifying currently used cropland was obtained from the spatial production allocation model (SPAM), whose spatial resolution was 10 km. The soil type data were obtained from the Harmonized World Soil Database (Nachtergaele et al., 2009).

Figure 1 Overall framework of emulator development for GGCMs. Each GGCM–crop combination was calibrated as an emulator independently. T: temperature, processed separately for daily maximum and minimum temperatures; P: precipitation; SRAD: solar radiation, N: nitrogen; soil: soil properties. When developing irrigated yield emulator, precipitation-related predictors are excluded.

2.2 Methods

Our study focused on the development and evaluation of emulators, which involves the following steps: (1) defining the predictors used to train the emulators, (2) preparing the predictors with climatic and non-climatic data, (3) training and cross-validating the emulators with a machine learning algorithm, and (4) evaluating the performance of emulators (Fig. 1).

Table 2 Predictors of emulation. For rainfed yield emulators, we used all these predictors, but for fully irrigated yield emulators, the precipitation predictors were not included. Full, GS, and MON were three strategies to develop emulators. Full: developing emulators with all the climatic predictors; GS: developing emulators with climatic predictors during the growing season; MON: developing emulators with climatic predictors on a monthly scale.

A checkmark in the row (√) indicates that the predictors were included in the emulator. The column “Time” clarifies the spatiotemporal dynamics of predictors. The number 1 represents both time- and space-variant predictors, 2 represents space-invariant predictors, and 3 represents time-invariant predictors.

2.2.1 Definition and preparation of predictors

All the predictors were computed or adapted from the GGCMs' input and output datasets. The climatic predictors were defined at two timescales, growing season (GS) and monthly (MON) (Table 2). The growing season average temperature, total precipitation, and average solar radiation were able to explain the variation of climatological mean yield of GGCMI phase 2 (Franke et al., 2020a). To improve the capacity of emulators to reproduce the year-to-year variation of crop model yield, daily variability and extremes of climate factors during the growing seasons were considered here. The variations of temperature, precipitation, and solar radiation during the growing seasons were calculated with the standard deviation of their daily values in each growing season, which represents the intensity of the daily fluctuation of weather. Additionally, heat and drought were selected to be the extreme climate predictors, which were quantified by extreme degree days (EDDs, cumulative temperature that exceed the high temperature threshold, Lobell et al., 2012) and maximum consecutive dry days (CDDs, maximum length of consecutive days without precipitation, Troy et al., 2015) because the negative effect of these two extremes can be shown by current GGCM (Heinicke et al., 2022). Other climate extremes, like excessive wetness, were not used because the GGCM failed to show the negative effect (Li et al., 2019a; Liu et al., 2022).

The monthly predictors only consisted of monthly average values. The monthly average temperature, total precipitation, and average solar radiation were harmonized according to the specific planting date. The number of months was determined with the crop-specific maximum growing season length over the global cropland defined by GGCMI phase 2. For winter and spring wheat, we prepared the climatic predictors over 10 and 8 months after sowing. For maize and rice, climatic predictors over 8 and 7 months after sowing were used, respectively.

The atmospheric CO₂ concentration and the nitrogen application rate were uniformly distributed predictors. All years and grid cells were set at the same CO₂ concentration and nitrogen application rate for each perturbation. Soil property is an important temporally constant predictor, whose interaction with climate played an important role in yield simulation and emulator development (Blanc, 2017). As the soil parameter settings of each GGCM varied, we selected the soil type at each grid to represent the spatial variation of soil properties. There were 13 soil types, including heavy clay, silty clay, light clay, silt clay loam, clay loam, silt, silt loam, sandy clay, loam, sandy clay loam, sandy loam, loamy sand, and sand. The most obvious difference across cultivars over global croplands is the growing degree requirement to reach maturity, which was determined by both mean climatology and cultivar traits. To reproduce the number of days from planting date to maturity date given by GGCMI phase 2 crop calendar input, we added a temporally constant growing season length as a predictor, i.e., a temporally constant growing season length.

As the purpose of emulator training is to develop a lightweight tool for crop simulation, there has always been a trade-off between the goodness of fit and the number of predictors. Therefore, we considered three strategies of using our predictors. Strategy A uses all predictors (the Full model), which is expected to derive the best goodness of fit. Strategy B uses only climatic predictors during the growing season (the GS model), together with CO₂ concentration, nitrogen application rate, site information, soil class, and growing season length. Strategy C uses only monthly average climatic predictors with other location-invariant predictors (the MON model). In general, Strategy B uses the smallest number of predictors, but those predictors need to be computed from daily climate forcing. Strategy C only relays on monthly climate data and is therefore the least costly strategy for data preparation. A comparison between the three strategies would help us find a good balance between the predictors used and overall goodness of fit of the emulators.

2.2.2 Emulator training and validation

Training the emulator of a specific GGCM is to derive the response relationship between input and output so that the emulator can replicate the complex process of yield simulation within the crop model. Emulation was trained by using machine learning regression on the GGCMI-2 ensemble of crop-specific simulated yield with all CTWN perturbations. Each grid year and perturbation combination was regarded as a sample in the fitting. We developed emulators of irrigated and rainfed yield and in A0 and A1 scenarios separately. Since the outputs of a GGCM outside current croplands were not well examined, we trained the machine-learning-based emulators only on currently used cropland, according to the SPAM data for each crop separately.

The extreme gradient boosting (XGBoost) algorithm, a highly efficient realization of the gradient boosting approach that showed the best performance in recent machine learning challenges (Chen and Guestrin, 2016), was used to train the emulators. Key parameters in XGBoost, including the learning rate (0.1), the number of estimators (4000), and the maximum tree depths (10), were tuned by a grid search along parameter dimensions based on the default parameter as a reference (Folberth et al., 2019). The goodness of fit of XGBoost was validated with the coefficient of determination $R_{\text{adjust}}^{\mathrm{2}}$:

$$R_{\text{adjust}}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{(n-\mathrm{1})\times (\mathrm{1}-R^{\mathrm{2}})}{n-k}},$$

where n is the sample size of the validation set, and k is the number of predictors.

Considering the spatiotemporal autocorrelation of simulated crop yield given by the GGCM, we now used a “held-out years and regions” strategy for leave-one-out cross-validation (Sweet et al., 2023). Specifically, all grid year samples are split into N folds. N is determined by the number of Köppen–Geiger (KG) classes, which have more than 100 grid cells with harvested areas. If there are too few harvested areas in one KG class, it will not be included in the cross-validation process. For each fold of emulator training and validation, we withhold 10 % of years (the last 3 years) and one entire KG class for validation, and the other grid year samples are used for training the emulator. We think selecting continuous years for validation can avoid temporal autocorrelation. If we randomly select 10 % of years, the correlation between adjacent years still exists. Actually, any continuous 3 years are able to solve this problem such that we just use the last years according to Sweet et al. (2023). Emulators were trained in Python3.8 with GPU (https://xgboost.readthedocs.io/en/latest/python/index.html, last access: 15 August 2023).

2.2.3 Evaluation of the emulator

Emulator performance was evaluated by comparing the 30-year emulated yield with the 30-year simulated yield of the GGCM. As we aimed at developing an emulator that could replicate the year-to-year variation of yield, the correlation coefficient (R), mean absolute error (MAE), and mean relative error (MRE) were used to evaluate the performance of emulators in the baseline and varied perturbations.

$$\begin{array}{c}{\displaystyle }R={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(Y_{\mathrm{XGB},i}-{\stackrel{\mathrm{‾}}{Y}}_{\mathrm{XGB}}\right)\left(Y_{\mathrm{GGCM},i}-{\stackrel{\mathrm{‾}}{Y}}_{\text{GGCM}}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left(Y_{\mathrm{XGB},i}-{\stackrel{\mathrm{‾}}{Y}}_{\mathrm{XGB}}\right)}^{\mathrm{2}}\cdot {\left(Y_{\mathrm{GGCM},i}-{\stackrel{\mathrm{‾}}{Y}}_{\mathrm{GGCM}}\right)}^{\mathrm{2}}}}}\\ {\displaystyle }\mathrm{MAE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\mathrm{|}{Y}_{\mathrm{XGB},i}-Y_{\mathrm{GGCM},i}\mathrm{|}}{n}}\\ {\displaystyle }\mathrm{MRE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\mathrm{|}(Y_{\mathrm{XGB},i}-Y_{\mathrm{GGCM},i})/Y_{\mathrm{GGCM},i}\mathrm{|}}{n}}\end{array}$$

where n is the sample size of the validation set, Y_GGCM,i is the annual simulated yield of the GGCMs, Y_XGB,i is the annual projected yield of the XGB algorithm, and ${\stackrel{\mathrm{‾}}{Y}}_{\mathrm{XGB}}$ and ${\stackrel{\mathrm{‾}}{Y}}_{\text{GGCM}}$ were the average XGBoost predicted and GGCM-simulated yield, respectively.

Table 3 Adjusted R² of XGBoost derived from 10-fold cross-validation with randomly selected samples.

The “–” symbol denotes no GGCM simulation; Winf: irrigated condition; W: rainfed condition. A0 denotes no adaptation, and A1 denotes cultivar adaptation to regain original growing season length under warming scenarios.

3 Results

3.1 Goodness of fit of the emulator training

Overall, the emulator developed with the XGBoost algorithm reproduced the variance of GGCM yield simulations well, with adjusted R² greater than 0.52 (Table 3). For most emulators the adjusted R² under a fully irrigated (Winf) simulation were greater than those under a rainfed simulation (W). Under A0 and A1 scenarios (A0 denotes no adaptation and A1 denotes adaptation of the growing season to regain the original growing season length under warming scenarios that otherwise lead to accelerated phenology and thus shorter growing seasons), the adjusted R² was comparable. For different crops, the performances of emulators developed for winter and spring wheat were slightly better than those developed for maize and rice. Among the GGCMs, PEPIC's behavior can best be emulated by emulators, with the greatest R² values for all crops and scenarios. There are also several GGCMs that are bit challenging for the XGB algorithm to capture, i.e., winter wheat and rice simulation from ORCHIDEE-crop, maize of pDSSAT, and spring wheat of LPJmL, with R² values ranging from 0.52 to 0.63.

Figure 2 Adjusted R² of emulators (10-fold cross-validation with randomly selected samples) with a different strategy used by the predictors. All: Full model, GS: GS model, MON: MON model. Emulators for ORCHIDEE by spring wheat and LPJ-GUESS by maize and rice were not fitted due to the lack of simulation by raw GGCMs.

The adjusted R² of emulators developed with all predictors (Full model) was greater than those developed with growing season predictors (GS model) and monthly predictors (MON model) (Fig. 2). GS models suffered from a reduced number of predictors, and their adjusted R² values were 0.1–0.15 smaller than corresponding MON models. Still, Full models had the largest adjusted R² at the cost of the largest number of predictors. For later usage of the emulators, a trade-off must be made between the cost of preparing predictors and model goodness of fit, and the MON model could be a balanced choice as it required only monthly average weather conditions.

Figure 3 Emulator performance in reproducing the year-to-year variation of global average yield (1981–2010) over current cropland. As ORCHIDEE-crop has not simulated yield under C360T0W0N200, we used the C360T0W10N200 as the baseline. Each point with the same color is yield in 30 years. R is the correlation coefficient and MAE is the mean absolute error.

3.2 Performance of emulators to capture the year-to-year variation of GGCM yield in the baseline

3.2.1 Performance of individual emulators at the global scale

Over current global cropland, the emulator of each GGCM could reproduce the year-to-year variation of global average yield well in the baseline period (during 1981–2010) (Fig. 3). All individual emulators captured the corresponding GGCM-simulated yield, with scatter concentrated on the 1:1 line. Different GGCM-simulated yield levels varied from 1.7 to 7.8 t ha⁻¹, but the performance of the emulators has not been influenced.

Figure 4 Correlation coefficient (R) and mean absolute error (MAE) over current cropland in the baseline (C360T0W0N200). As ORCHIDEE-crop has not simulated yield under C360T0W0N200 perturbation, we used the C360T0W10N200 as the baseline. The dot denotes the median and the error bar denotes the interquartile range from all grid cells in which the crop is grown according to SPAM 2010 data.

3.2.2 Performance of individual emulators at the grid scale

The overall performances of emulators at grid level were good for most crop–GGCM combinations in the baseline. The performances of each emulator over current global cropland grids were plotted by using the scatter of MAE and R (Fig. 4). The capacity of emulators to reproduce the wheat yield simulated by GGCMs was better than that of maize and rice. The median R values over current winter and spring wheat harvested areas were greater than 0.7. The R of the EPIC-TAMU emulator and the LPJ-GUESS emulator were relatively smaller than the other eight emulators developed for winter and spring wheat, respectively. The median MAEs over current winter and spring wheat harvested areas were less than 0.4 and 0.3 t ha⁻¹ for winter and spring wheat emulators, respectively, and the MAEs of the pDSSAT emulator and LPJmL emulator were relatively greater. Over current maize harvested areas, the median R was greater than 0.6 and the median of MAE was less than 0.7 t ha⁻¹, except the pDSSAT emulator. The median R values of emulators developed for rice were greater than 0.5, and the median MAEs were less than 0.4 t ha⁻¹ over current rice harvested areas, whereas the performances of the pDSSAT emulator and CARAIB emulator were relatively worse.

Figure 5 Multi-model ensemble median R in the baseline over current cropland. R: correlation coefficient between simulated and emulated yield time series of each GGCM from 1981 to 2010.

3.2.3 Performance of multiple-emulator ensemble at the grid scale

The multi-emulator ensemble median was able to reproduce the year-to-year variation of gridded yield over current cropland in the baseline (C360T0W0N200) from 1981 to 2010. The temporal correlation coefficients R between GGCM-simulated and emulated yield time series over most current harvested areas were greater than 0.7 (multi-model ensemble median) (Fig. 5), and the uncertainty (standard deviation) of R across emulators was smaller than 0.3 (Fig. S1 in the Supplement). The mean absolute error (MAE) and mean relative error (MRE) over most current harvested areas were mostly smaller than 1 t ha⁻¹ and 30 %, respectively (Fig. S2 in the Supplement). The spatial pattern of MRE for four crops all showed a hotspot of large MRE in the Middle East, and for maize the hotspot of large MRE was also found in southern China (Fig. S2).

Figure 6 Performance of one exemplary emulator (LPJmL-A0) in reproducing the year-to-year variation of global mean yield from 1981 to 2010 under varied individual CTWN perturbations. Each point with the same color is yield in 1 year. The performances of other emulators are similar to LPJmL-A0, which can be seen Tables S1 and S2.

3.3 Performance of emulators to capture the year-to-year variation of GGCM yield in the CTWN cube

3.3.1 Performance of individual emulators at the global scale

The agreement of year-to-year variation of global average yield between simulation and emulation was consistent with changes in CTWN cubes over present cropland (Fig. 6). Under varied CTWN perturbations, the emulator could reproduce the year-to-year variation of global mean yield from 1981 to 2010 well. Even when the temperature perturbation reached +6 K, the emulator was still able to capture the year-to-year variation of global mean yield. Similarly, when the precipitation was less than the baseline by 50 %, the year-to-year variation of emulated global mean yield was well matched with those of GGCM simulations. Additionally, fertilization with elevated CO₂ concentrations and nitrogen application have been well reproduced by the emulator. Similar capacity to reproduce the annual global mean yield was also found in other emulators (Tables S1 and S2 in the Supplement). Even under concurrent warm and drought conditions, i.e., T + 6 K and W − 50 %, the year-to-year variation of global mean yield could be well reproduced by the emulator (Fig. S3 in the Supplement).

Figure 7 R–MAE scatter of the exemplary emulator (LPJmL-A0) under varied single CTWN perturbations. Each dot denotes the median of R or MAE over current cropland, and the error bar denotes the interquartile range. R: correlation coefficient, MAE: mean absolute error. More details of other emulators can be seen in Tables S1 and S2.

3.3.2 Performance of individual emulators at the grid scale under single perturbation

To illustrate the performance of individual emulators in reproducing annual yield variation, we selected the LPJmL-A0 emulator as an example. The R–MAE scatter plots of LPJmL-A0 illustrate the response of gridded accuracy to varied perturbations of CTWN (Fig. 7). The changes in the accuracy of emulators under single CTWN perturbations were small, with the largest differences in spring wheat for modifications in the CO₂ (C) and nitrogen (N) dimensions. The overall accuracy could be kept at a high level, with greater R and smaller MAE. Under temperature perturbations, the median R values of emulators for four crops were greater than 0.7, and the range of R values was smaller than 0.2. The median MAEs of emulators for four crops were less than 0.55, and the variation of median MAEs was smaller than 0.2 from +1 to +6 K perturbations. For precipitation perturbations, the median R values of emulators for four crops were greater than 0.85; meanwhile, the difference of median R values across varied precipitation perturbations was smaller than 0.1. The median MAEs of emulators for four crops was smaller than 0.3, and the range of median MAEs variation was as small as 0.06. The median R values and MAEs of emulators for four crops under CO₂ concentration perturbations and nitrogen perturbations were comparable to those under temperature and precipitation perturbations, except for spring wheat. Although the performance of the spring wheat emulator under CO₂ and nitrogen perturbations was not as good as other crops, the median R values were still greater than 0.75, and the median MAEs were smaller than 0.6. Similar patterns of other emulator performances under single perturbations at the grid scale are shown in Tables S1 and S2.

Figure 8 Correlation coefficient (R) of the multiple-emulator ensemble under varied TW perturbations. The line denotes the median of R over current cropland, and the shaded area denotes the range of median R over current cropland across emulators.

Figure 9 Mean absolute error (MAE) of the multiple-emulator ensemble under varied CTWN perturbations. The line denotes the median of R over current cropland, and the shaded area denotes the range of median R over current cropland across emulators.

3.3.3 Performance of multiple-emulator ensemble at the grid scale under single perturbation

When looking at the ensemble of multiple emulators, the Rs and MAEs under CTWN cubes was not divergent obviously (Figs. 8 and 9).

Under temperature perturbations, the range of model ensemble median Rs across multiple emulators was smaller than 0.2, and the range of median MAEs was as small as 0.4 t ha⁻¹. For precipitation perturbation, the difference in median Rs was less than 0.03, and the changes in median MAEs were less than 0.1 t ha⁻¹. Under the perturbation of CO₂ concentration, the emulators for winter wheat, maize, and rice showed greater median R values which ranged from 0.89 to 0.98. The variation of median MAEs was smaller than 0.09 t ha⁻¹. The median R values of the emulator for spring wheat, however, tended to substantially decline under 810 ppm perturbation, and the median MAEs tended to increase simultaneously. Similarly, for nitrogen perturbation, the range of median R values was less than 0.27, and the range of median MAEs was smaller than 0.3 t ha⁻¹, except for emulators of spring wheat and rice. The declined R and increased MAE were caused by the reduction of valid sample size from the GGCM output yield under CO₂ and nitrogen perturbations (Figs. S4 and S5 in the Supplement).

Figure 10 Fraction of harvested areas with a high correlation coefficient (R > 0.8) under varied T + W perturbations. Example of EPIC-TAMU-A0, LPJmL-A0, and pDSSAT-A0 emulators because only these three GGCMs contain the full range of CTWN perturbations for all four crops.

Figure 11 Fraction of harvested areas with low mean absolute error (MAE < 0.5 t ha⁻¹) under varied T + W perturbations. Example of EPIC-TAMU-A0, LPJmL-A0, and pDSSAT-A0 emulators because these three GGCMs contain the full range of CTWN perturbations for all four crops.

3.3.4 Performance of multiple emulators at grid the scale under dual perturbations

The performance of emulators was influenced by changes in simultaneous perturbations in two different CTWN dimensions (dual perturbations). The emulators performed well over most current cropland, but at extreme increases in T and reductions in W (Fig. 10), the emulators could represent the GGCMI-simulated year-to-year variation only on substantially smaller shares of the current cropland. The fraction of current areas with R greater than 0.8 was the highest in the baseline but decreases under warmer and drier conditions. The fraction was reduced to less than 40 % under the compound T +6 K and W −50 % perturbation, which illustrates the poor capacity of the emulator under compound hot–dry conditions. However, the fraction of harvested areas with MAE smaller than 0.5 t ha⁻¹ did not vary much across T + W perturbations (Fig. 11). The performances of emulators under dual perturbations for wheat were better than those for maize and rice. The fraction of maize and rice harvested area with R greater than 0.8 was relatively smaller than that of wheat. The maize harvested area with MAE smaller than 0.5 t ha⁻¹ was smaller than other crops. Among the three GGCMs with the full range of CTWN perturbations, the fraction of harvested area with high accuracy for the LPJmL emulator and the pDSSAT emulator was more than EPIC-TAMU emulator.

4 Discussion

4.1 Emulator trained to capture year-to-year variation in crop yield

Our emulator was designed to reproduce the year-to-year variation of crop yield. Therefore, the annual yield was the target variable in emulator fitting. To capture the year-to-year crop yield variation well, the climatic predictors were divided into growing season average, daily variation, and climatic extremes to capture the possible drivers of yield variation. Predictor engineering referred to the existing knowledge compiled in crop models that year-to-year variation of crop yield is associated with growing season temperature and precipitation (Ray et al., 2015), extreme heat (Iizumi and Ramankutty, 2016), and drought (Heinicke et al., 2022). Temperature and precipitation have been confirmed to be the dominant drivers of crop yield variability (Schauberger et al., 2016). Moreover, the interaction between soil type and climate was considered in our emulator design. Although CO₂ concentration and soil type were not regarded as important contributors to yield variability, their interaction with climate could also influence the yield variability (Kadam et al., 2014). The role of soil type has been revealed by previous emulators fitted by multivariate regression; the average effect of temperature and precipitation differed greatly depending on soil type (Blanc, 2017). Compared with the emulator designed to reproduce the climatological mean yield, our emulator is more suitable to project changes in yield variability (Liu et al., 2021b).

We developed the emulators with one statistical relationship for each crop between GGCM-simulated yield and predictors for all grids over global land. Each grid cell represents a sample in the soil–climate–fertilizer continuum, and the training data have no lateral relationships. However, the response of simulated crop yield to climate change was spatially heterogeneous, which mainly depends on the cultivars. Therefore, one statistical relationship between yield and climatic predictors was not fully appropriate for each grid. In response, we used the length of growing season, a representative predictor of cultivar characteristics, to adjust the global statistical relationship to each grid. Therefore, predictors contained both temporal varied and constant variables. The temporally varied predictors were climatic variables which mainly played a role in reproducing the annual yield variation, and the temporally constant predictors were non-climatic variables, like growing season length, that delineated the spatial distinction of crop yield response to climate. Compared with region-specific emulator development, combining the temporally varied and constant predictors was more concise and could profit from a broader range of data in the training set.

4.2 Potential application of the emulators performing well in related fields

The good performance over most grid cells indicated the potential capacity of emulators in spatiotemporal downscaling, projecting annual yield in the future, and multi-model ensemble simulation.

The emulator could be used to conduct spatiotemporal yield downscaling because of the good performance of individual emulators in reproducing the annual crop yield variation in the baseline. As the emulator in this study was developed with a regression-based machine learning technique by using all the grid year data points, the emulation is not limited to the spatial resolution of the training data. The emulator can be applied to any point with information on the predictors and can produce yield projections as finely resolved as the forcing input. From the aspect of time series of yield, the raw GGCM data include empty values (NaN – not a number) in some year and grid cell data points, which may be caused by the lack of regional data for calibration. The vacancy of yield time series in some grids could be imputed by the emulator (Folberth et al., 2019), similar to studies which generated spatiotemporally continuous gridded crop yield data (Chen et al., 2022; Iizumi et al., 2014).

The emulator was able to project the annual yield in the future climate scenarios, which depends on the individual emulator performing well in reproducing annual yield under CTWN cubes. In contrast to many previous emulators developed with historical crop model simulations (Xu et al., 2021), our emulator could reproduce the CO₂ fertilization effect, which is an important forcing for the future. The recently developed emulator based on GGCMI phase 2 simulation under CTWN cubes could only project the climatological mean yield because the target variable in the emulation was the climatological mean yield (Franke et al., 2020a). In contrast, our emulator can project the annual yield variation and is not constrained by the maximum warming considered in the GGCMI phase 2 dataset (T +6 K), but by the maximum temperature within the training dataset (warmest grid cell +6 K), so that the applicability is broader (Müller et al., 2021).

It is more efficient to conduct multi-model ensemble simulations with emulators than GGCMs, as the emulators show good skill in reproducing GGCM results, and the emulators drastically reduce the computational time, memory requirement, and expertise to operate process-based crop models. First, the input of multiple emulators was consistent and compatible, but the inputs of raw GGCMs were divergent and incompatible because the structure of input data and file format of each GGCM were designed independently. Second, the timescale of emulator input was monthly or growing seasonal, which was less complex than daily inputs of GGCMs. Apart from the ensemble simulation, the multiple emulators could also be used to explore and disentangle the uncertainty across models.

4.3 Uncertainties

The weaknesses of machine learning algorithms and raw GGCMs have brought some uncertainties into the emulators. The uncertainties induced by the machine learning algorithm were as follows.

1. When the climate factors went beyond the range of training data, the weakness of machine learning in out-of-sample prediction could bring great uncertainty. The emulator inputs should be capped by the range of training data. The limit of our emulator was the warmest grid under +6 K perturbation. As there is polar amplification, the strongest warming mostly happens in cooler regions. Thus, the projected temperature exceeding the training range would not be widespread over global croplands.

2. Although the emulators could reproduce the GGCM annual yield with high accuracy in most cases, there were cases in which the machine learning algorithm did not show good reproduction skill. As the emulated function was intended to smooth the response of simulated crop yield to climate, samples at the margins of training data tend to show lower emulator skill. Extreme conditions, i.e., +6 K, −50 % water, 810 ppm, 10 kg ha⁻¹, show reduced R and increased MAE. Using the emulators to estimate annual crop yield under extreme perturbation conditions should be done with caution and the additional uncertainty induced by the emulators should be considered in the interpretation of results.

3. Last but not least, as the emulators are intended to be lightweight tools that can replicate raw GGCMs, their capability to simulate crop yields is limited to the capability of the original GGCMs. This raises the issue that emulators are unlikely to show good performance in simulating crop yield responses to climate extremes, exactly like the raw GGCMs, which have shown poor performance in capturing the yield impact of heatwave and drought (Heinicke et al., 2022), and the lack of negative effects of excessive wetness (Li et al., 2019a). Resolving such a problem requires the improvement of raw GGCMs' capability to simulate yield response to climate extremes or statistical promotion of the GGCMs' outputs under extreme weather events.

5 Conclusion

In this study, we developed machine-learning-based statistical crop yield emulators to reproduce the year-to-year variation of crop yield in response to perturbations in CO₂ concentration, temperature, water, and nitrogen application rate from the GGCMI phase 2 archives. To examine the potential value of these emulators, we evaluated the performance of emulators at the global and gridded scale under a baseline with single and dual perturbations.

The results indicated that the performance of emulators was good enough to reproduce the year-to-year variation of global average crop yield in the baseline (R > 0.9), and the difference in accuracy between individual GGCM emulators was not large. Similarly, under single and dual perturbations, the capacity of emulators to reproduce the year-to-year variation of global mean crop yield was not substantially changed. At the gridded level, the performance of emulators over most of the current croplands in the baseline was still good in the sense that R was greater than 0.6 and MAE was smaller than 1 t ha⁻¹. The performance of individual emulators was consistently good under single CTWN perturbations, without substantial changes in R and MAE. Similarly, multiple emulators also performed well in reproducing the annual yield under single CTWN perturbations, and most grid cells across the current cropland showed greater R and smaller MAE under simultaneous perturbations of T and W. The overall good capacity of emulators to reproduce the year-to-year variation of GGCM-simulated crop yield indicated the role of emulators in spatiotemporal downscaling, crop yield projection, and multi-model ensemble simulation. The emulators were able to boost the ability to assess crop yield failure risk under future climate change and help to better understand food stability and climate risk adaptation.