Simulation of crop yield using the global hydrological model H08 (crp.v1)

A better understanding of the food–water nexus requires the development of an integrated model that can simultaneously simulate food production and the requirements and availability of water resources. H08 is a global hydrological model that considers human water use and management (e.g., reservoir operation and crop irrigation). Although a crop growth sub-model has been included in H08 to estimate the global crop-specific calendar, its performance as a yield simulator is poor, mainly because a globally uniform parameter set was used for each crop type. In addition, the effects of CO₂ fertilization and vapor pressure deficit on crop yield were not considered. Here, through country-wise parameter calibration and algorithm improvement, we enhanced H08 to simulate the yields of four major staple crops: maize, wheat, rice, and soybean. The simulated crop yield was compared with the Food and Agriculture Organization (FAO) national yield statistics and the global dataset of historical yield for major crops (GDHY) gridded yield estimates with respect to mean bias (across nations) and time series correlation (for individual nations). Our results showed that the effects of CO₂ fertilization and vapor pressure deficit had opposite impacts on crop yield. The simulated yield showed good consistency with FAO national yield. The mean biases of the major producer countries were considerably reduced to 2 %, 2 %, −2 %, and −1 % for maize, wheat, rice, and soybean, respectively. The capacity of our model to capture the interannual yield variability observed in FAO yield was limited, although the performance of our model was comparable to that of other mainstream global crop models. The grid-level analysis showed that our model showed a similar spatial pattern to that of the GDHY yield in terms of reproducing the temporal variation over a wide area, although substantial differences were observed in other places. Using the enhanced model, we quantified the contributions of irrigation to global food production and compared our results to an earlier study. Overall, our improvements enabled H08 to estimate crop production and hydrology in a single framework, which will be beneficial for global food–water nexus studies in relation to climate change.

1 Introduction

Food security has become an important global challenge because of the growing population and increasing competition for crop usage (Ray et al., 2022). A key factor in food security is crop production, which is largely affected by irrigation water availability, particularly in regions with insufficient precipitation (Chiarelli et al., 2022). Currently, for example, approximately 40 % of global crop production relies on irrigation (Perrone, 2020). The use of water for this irrigation causes approximately 65 % of global total water withdrawal and 90 % of global water consumption (Shiklomanov, 2000; Döll and Siebert, 2002). These high rates of withdrawal and consumption have negative consequences for both surface water and groundwater systems, such as river fragmentation and groundwater table decline (McDermid et al., 2021; Perrone, 2020). To minimize such negative consequences, there is an increasing impetus toward sustainable water use (McDermid et al., 2021; Perrone, 2020; Rosa et al., 2018, 2020; Okada et al., 2018; Ai et al., 2021). To more fully address the complex interactions between crop production and sustainable water management, accurate representations of crop growth and water cycle with human activities should be placed within a consistent model framework during the development of an integrated model.

Many models have successfully incorporated the crop growth process and can simulate the global crop yield. These include LPJmL (Bondeau et al., 2007; Fader et al., 2010), GEPIC (Liu et al., 2007), PEGASUS (Deryng et al., 2011), CLM-Crop (Drewniak et al., 2013), PRYSBI2 (Sakurai et al., 2014), pAPSIM (Elliott et al., 2014), pDSSAT (Elliott et al., 2014), CROVER (Okada et al., 2015), ORCHIDEE-crop (Wu et al., 2016), PEPIC (Liu et al., 2016), MATCRO (Masutomi et al., 2016), SIMPLACE-LINTUL5 (Webber et al., 2016), PROMET (Zabel et al., 2019; Degife et al., 2021), and ACEA (Mialyk et al., 2022). However, only a few of these models, such as LPJmL and CROVER, have globally implemented schemes for irrigation constrained by detailed spatiotemporal water availability (i.e., explicit consideration of river routing and water withdrawal). The lack of inclusion of such schemes severely limits the ability of these models to be used in comprehensive investigations of global food–water tradeoffs, particularly in terms of specifying the sources of water withdrawal used for crop irrigation.

In this study, we developed a new crop–water global model based on the H08 global hydrological model (Hanasaki et al., 2008a, 2018). Although H08 has detailed functions for specifying water sources and estimating crop specific yield based on the formulations of the SWAT model (Neitsch et al., 2002), its performance as a crop yield simulator has been poor in comparison with the Food and Agriculture Organization (FAO) yield statistics and other gridded yield datasets. This poor performance is mainly because of the adoption of the global uniform parameters related to crop growth. These default parameters are acquired from the SWIM model, a variant of the SWAT model (Arnold et al., 1994), which is mainly for use in Europe and temperate climate zones (Krysanova et al., 2000). This leads to overestimation or underestimation when it is used in other regions with different crop management practices and climatic conditions. Additionally, the effects of CO₂ fertilization (Stockle et al., 1992) and changes in vapor pressure deficit (Stockle and Kiniry, 1990) on crop yield have not yet been considered. These two factors are particularly important when analyzing the impacts of climate change on crop yield (Jägermeyr et al., 2021; Yuan et al., 2019).

Despite multiple attempts to optimize the parameters involved, global crop yield simulation remains challenging. For example, Fader et al. (2010) proposed the concept of management intensity, which represents the degree and frequency of field agronomy management (e.g., fertilizer, technology, and weed control). They adopted this concept in a global vegetation model, LPJml, by adjusting a key parameter of maximum leaf area index at the country level, which exhibited good agreement between the calibrated yield and FAO yield statistics. This adjustment enabled LPJml to be used in investigations of the crop–water relations by estimating crop water productivity and virtual water content (Fader et al., 2010). Deryng et al. (2011) calibrated the light use efficiency coefficient based on spatially explicit crop yield data reported by Monfreda et al. (2008). Iizumi et al. (2009) developed a large-scale crop model for paddy rice in Japan, known as the PRYSBI model, whereby multiple parameters were calibrated via the Markov chain Monte Carlo technique at the subnational level. The results showed that the Markov chain Monte Carlo method is a powerful approach for optimizing multiple parameters in a nonlinear and complex model. Sakurai et al. (2014) used a similar method globally and estimated eight parameters based on Free-Air Carbon Dioxide Enrichment (FACE) data with hundreds of thousands of calculation steps in the Markov chain Monte Carlo process. Each of the above methods has its own advantages and disadvantages. For example, the method of Fader et al. (2010) was based on FAO national yield statistics, whereas the methods in the other three studies require spatially explicit yield data. Additionally, Fader et al. (2010) and Deryng et al. (2011) mainly focused on a single parameter, whereas Iizumi et al. (2009) and Sakurai et al. (2014) addressed multiple parameters.

To enhance the capacity of H08 to simulate the yields of four major staple crops (i.e., maize, wheat, rice, and soybean), we first added two new functions to the H08 crop sub-model by considering the effects of CO₂ fertilization and vapor pressure deficit change on crop yield. Then, we adopted the method of Fader et al. (2010) for parameter calibration because of its robust performance, minimal computation costs, and simplicity of implementation, as well as because the method requires only national yield data which are easily accessible and generally reliable. Next, we evaluated model performance with respect to mean bias, time series variation, and time series correlation in accordance with the general framework proposed by Müller et al. (2017) using FAO statistical national data and recently published grid-level data. We sought to determine how crop yield responds to changes in CO₂ concentration and vapor pressure deficit to determine whether the enhanced H08 model could reproduce the mean historical yield at the national scale and to determine whether the model could also capture interannual variation in historical yield times series; we also aimed to compare spatial time series correlations with other spatially explicit data. Finally, we investigated the contributions of irrigation to the global production of maize, wheat, rice, and soybean using the enhanced model as a case study for its application.

2 Materials and methods

2.1 H08 overview

H08 is a global hydrological model that includes natural and anthropogenic hydrological processes at a spatial resolution of 0.5^∘ and a temporal resolution of 1 d. It was developed with six sub-models: land surface hydrology, river routing, crop growth, reservoir operation, environmental flow requirements, and anthropogenic water withdrawal (Hanasaki et al., 2008a). It has been updated with several new schemes including groundwater recharge and abstraction, aqueduct water transfer, local reservoirs, seawater desalination, and return flow and delivery loss (Hanasaki et al., 2018). With these newly added functions, H08 is one of the most detailed global hydrological models available for the estimation of sector-wise and water-source-wise water withdrawal and availability. In the agriculture sector, H08 can estimate irrigation water demand and supply on a daily and grid-cell basis with several unique features. First, it can estimate the irrigation water withdrawal from both renewable and non-renewable groundwater sources. Second, it considers the effects of irrigation water withdrawal in the upper stream. Third, it includes the influence of reservoir operation on irrigation water availability. H08 was fully described in multiple studies (Hanasaki et al., 2008a, b, 2018).

Figure 1 Comparison of the mean yield from 1986 to 2015 of different simulations and FAO statistics. (a) Maize, (b) wheat, (c) rice, and (d) soybean. Further details on the five utilized simulations (D, C, V, CV, and CVC) are listed in Table 1.

2.2 Crop sub-model

2.2.1 Overview

The crop growth sub-model accumulates plant biomass at a daily interval until physiological maturity; it also simulates phenological development. The daily increase in potential biomass (ΔB) (kg ha⁻¹) is estimated based on radiation use efficiency and photosynthetic active radiation using the method of Monteith et al. (1977) (Eq. 1). Crop phenological development is based on daily heat unit accumulation theory, whereby physiological maturity is reached when the accumulated daily heat unit value is equal to the potential heat unit value. The harvest index is used to partition the total aboveground biomass with respect to grain yield. Regulating factors, including water and air temperature, are used to adjust the yield variation. A schematic figure that shows the basic biophysical processes of the crop sub-model is shown in Fig. 1b in Ai et al. (2020). Although the algorithm is based on SWAT and SWIM, and a detailed description was previously provided (Hanasaki et al., 2008a; Ai et al., 2020), the main formulation is briefly described below because it is an important foundation for the forthcoming discussion on parameter calibration.

2.2.2 Basic algorithms

ΔB is calculated as follows:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }B=\mathrm{be}\cdot \mathrm{PAR}\cdot \mathrm{REGF},\end{array}$$

where be is a crop-specific parameter of radiation use efficiency, PAR is photosynthetically active radiation, and REGF is the crop regulating factor. PAR is calculated using shortwave radiation (R_s) (W m⁻²) and leaf area index (LAI) as follows:

$$\begin{array}{}\text{(2)}& \mathrm{PAR}=\mathrm{0.02092}\cdot R_{\mathrm{s}}\cdot [\mathrm{1}-\exp (-\mathrm{0.65}\cdot \mathrm{LAI}\left)\right].\end{array}$$

LAI is calculated according to the growth stage indicated by the heat unit index (I_hun), which is calculated as the ratio of accumulated daily heat units ∑H_una(t) and the potential heat unit (H_un):

$$\begin{array}{}\text{(3)}& I_{\text{hun}}={\displaystyle \frac{\sum H_{\text{una}}\left(t\right)}{H_{\text{un}}}}.\end{array}$$

The daily heat units H_una(t) are based on the difference between the daily mean air temperature (T_a) and the crop's specific base temperature (T_b; provided as a crop-specific parameter):

$$\begin{array}{}\text{(4)}& H_{\text{una}}\left(t\right)=T_{\mathrm{a}}-T_{\mathrm{b}},\end{array}$$

if $I_{\text{hun}}<\lfloor \mathrm{dpl}\mathrm{1}\rfloor \cdot \mathrm{0.01}$,

$$\begin{array}{}\text{(5)}& \mathrm{LAI}={\displaystyle \frac{(\mathrm{dpl}\mathrm{1}-\lfloor \mathrm{dpl}\mathrm{1}\rfloor )\cdot I_{\text{hun}}}{\lfloor \mathrm{dpl}\mathrm{1}\rfloor \cdot \mathrm{0.01}}}\cdot \mathrm{blai},\end{array}$$

if $\lfloor \mathrm{dpl}\mathrm{1}\rfloor \cdot \mathrm{0.01}\le I_{\text{hun}}<\lfloor \mathrm{dpl}\mathrm{2}\rfloor \cdot \mathrm{0.01}$,

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{LAI}=& \mathit{\left\{}\right(\mathrm{dpl}\mathrm{1}-\lfloor \mathrm{dpl}\mathrm{1}\rfloor )\\ & +{\displaystyle \frac{\begin{array}{c}\left[(\mathrm{dpl}\mathrm{2}-\lfloor \mathrm{dpl}\mathrm{2}\rfloor )-(\mathrm{dpl}\mathrm{1}-\lfloor \mathrm{dpl}\mathrm{1}\rfloor )\right]\\ \times (I_{\text{hun}}-\lfloor \mathrm{dpl}\mathrm{1}\rfloor \cdot \mathrm{0.01})\end{array}}{\lfloor \mathrm{dpl}\mathrm{2}\rfloor \cdot \mathrm{0.01}-\lfloor \mathrm{dpl}\mathrm{1}\rfloor \cdot \mathrm{0.01}}}\mathit{\}}\cdot \mathrm{blai},\end{array}\end{array}$$

if $\lfloor \mathrm{dpl}\mathrm{2}\rfloor \cdot \mathrm{0.01}\le I_{\text{hun}}<\mathrm{dlai}$,

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{LAI}=& \mathit{\left\{}\right(\mathrm{dpl}\mathrm{2}-\lfloor \mathrm{dpl}\mathrm{2}\rfloor )\\ & +{\displaystyle \frac{\begin{array}{c}\left[\mathrm{1}-(\mathrm{dpl}\mathrm{2}-\lfloor \mathrm{dpl}\mathrm{2}\rfloor )\right]\\ \times (I_{\text{hun}}-\lfloor \mathrm{dpl}\mathrm{2}\rfloor \cdot \mathrm{0.01})\end{array}}{\mathrm{dlai}-\lfloor \mathrm{dpl}\mathrm{2}\rfloor \cdot \mathrm{0.01}}}\mathit{\}}\cdot \mathrm{blai},\end{array}\end{array}$$

if dlai<I_hun,

$$\begin{array}{}\text{(8)}& \mathrm{LAI}=\mathrm{16}\cdot \mathrm{blai}(\mathrm{1}-I_{\text{hun}}{)}^{\mathrm{2}},\end{array}$$

where dlai is the fraction of growing season when growth declines, dpl1 and dpl2 are shape parameters of the LAI growth curve (see the definition in Table 1 in Ai et al., 2020), and blai is the maximum leaf area index.

Table 1 Simulation settings.

REGF is calculated as

$$\begin{array}{}\text{(9)}& \mathrm{REGF}=min(T_{\mathrm{s}},W_{\mathrm{s}},N_{\mathrm{s}},P_{\mathrm{s}}),\end{array}$$

where T_s, W_s, N_s, and P_s are the stress factors for temperature, water, nitrogen, and phosphorous, respectively. The details of water and temperature stress are provided in the work of Ai et al. (2020). Nitrogen and phosphorous stress were not considered in the original model (Hanasaki et al., 2008a) and were indirectly represented in the calibration simulation in the present study.

The aboveground biomass (B_ag) (kg ha⁻¹) is estimated with the accumulated biomass (∑ΔB) as

$$\begin{array}{}\text{(10)}& B_{\text{ag}}=[\mathrm{1}-(\mathrm{0.4}-\mathrm{0.2}\cdot I_{\text{hun}}\left)\right]\sum \mathrm{\Delta }B.\end{array}$$

The crop yield (Yld) (kg ha⁻¹) is finally estimated from the aboveground biomass (B_ag) using the crop-specific harvest index (Harvest) on the date of the harvest:

$$\begin{array}{}\text{(11)}& \mathrm{Yld}=\mathrm{Harvest}\cdot {\displaystyle \frac{\mathrm{WSF}}{\mathrm{WSF}+\exp (\mathrm{6.117}-\mathrm{0.086}\cdot \mathrm{WSF})}}\cdot B_{\text{ag}},\end{array}$$

where WSF is the ratio of SWU (accumulated actual plant evapotranspiration in the second half of the growing season) to SWP (accumulated potential evapotranspiration in the second half of the growing season):

$$\begin{array}{}\text{(12)}& \mathrm{WSF}={\displaystyle \frac{\mathrm{SWU}}{\mathrm{SWP}}}\cdot \mathrm{100}.\end{array}$$

Differences in crop type are expressed by the differences in crop parameters (e.g., be, blai, and T_b). Currently, the crop sub-model can simulate the yield for 18 food crops. The globally uniform default parameters for the food crops were collected from the default parameters of the SWIM model (Krysanova et al., 2000).

2.3 Algorithm improvement

Here, the crop sub-model was improved as follows. First, the effects of CO₂ fertilization and vapor pressure deficit change on radiation use efficiency were added to the H08 crop sub-model using the equations and parameters adopted in SWAT (Neitsch et al., 2011; Arnold et al., 2013). Specifically, the radiation use efficiency (be) is adjusted according to the concentration of CO₂ as

$$\begin{array}{}\text{(13)}& \mathrm{be}={\displaystyle \frac{\mathrm{100}\cdot {\mathrm{CO}}_{\mathrm{2}}}{{\mathrm{CO}}_{\mathrm{2}}+\exp (r_{\mathrm{1}}-r_{\mathrm{2}}\cdot {\mathrm{CO}}_{\mathrm{2}})}},\end{array}$$

where be is the radiation use efficiency, CO₂ is the CO₂ concentration in the atmosphere (ppmv), and r₁ and r₂ are shape coefficients defined as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }{r}_{\mathrm{1}}=\ln ⌊{\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}\mathrm{amb}}}{\mathrm{0.01}\cdot {\mathrm{be}}_{\mathrm{amb}}}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{amb}}⌋+r_{\mathrm{2}}\cdot {\mathrm{CO}}_{\mathrm{2}\mathrm{amb}},\text{(15)}& {\displaystyle }{r}_{\mathrm{2}}={\displaystyle \frac{\ln ⌊\frac{{\mathrm{CO}}_{\mathrm{2}\mathrm{amb}}}{\mathrm{0.01}\ast {\mathrm{be}}_{\mathrm{amb}}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{amb}}⌋-\ln ⌊\frac{{\mathrm{CO}}_{\mathrm{2}\mathrm{hi}}}{\mathrm{0.01}\cdot {\mathrm{be}}_{\mathrm{hi}}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{hi}}⌋}{{\mathrm{CO}}_{\mathrm{2}\mathrm{hi}}-{\mathrm{CO}}_{\mathrm{2}\mathrm{amb}}}},\end{array}$$

where CO_2amb is the ambient atmospheric CO₂ concentration (ppmv), CO_2hi is an elevated atmospheric CO₂ concentration (ppmv), be_amb is the be of the crop at CO_2amb, and be_hi is the be of the crop at CO_2hi.

Additionally, the be is adjusted with the vapor pressure deficit (vpd) (kPa) as

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}\mathrm{be}=& {\mathrm{be}}_{\mathrm{vpd}=\mathrm{1}}-\mathrm{\Delta }{\mathrm{be}}_{\mathrm{dcl}}\\ & \times (\mathrm{vpd}-{\mathrm{vpd}}_{\mathrm{thr}})\text{if }\mathrm{vpd}>{\mathrm{vpd}}_{\mathrm{thr}},\end{array}\text{(17)}& {\displaystyle }\mathrm{be}={\mathrm{be}}_{\mathrm{vpd}=\mathrm{1}}\text{if }\mathrm{vpd}\le {\mathrm{vpd}}_{\mathrm{thr}},\end{array}$$

where be_vpd=1 is the be for the plant at a vpd of 1 kPa, Δbe_dcl is the rate of be decline per unit increase in vpd, and vpd_thr is the threshold vpd above which a plant will exhibit reduced radiation use efficiency. The vpd_thr value is assumed to be 1 kPa.

2.4 Parameter calibration

Next, we calibrated the key parameter of maximum leaf area index (blai) and adjusted the harvest index (Harvest) accordingly by adopting the concept of management intensity in accordance with the method of Fader et al. (2010). For many countries in the world, the historical annual crop yield from FAO data shows an increasing trend. Hence, the common method of splitting data into two periods, one for calibration and one for validation, was not applicable here. Therefore, we used the mean of even years for calibration and the mean of odd years for confirmation. Specifically, we calibrated the maximum leaf area index by iterating the values from 0.5 to 7.1, with an interval of 0.3, under both rainfed and irrigation conditions in the even years from 1986 to 2015. The crop-specific best maximum leaf area index in each country was then determined as the value that can minimize the bias between the mean simulated yield and mean FAO statistical yield. When FAO statistical yield or simulated yield data were missing for a country, we used the original crop-specific default values. Then, we adjusted the harvest index with the calibrated blai (see Table S1 in the Supplement). The calibration and confirmation results showed good agreement with the FAO statistics (Fig. S1 in the Supplement).

2.5 Meteorological data

The ISIMIP3a GSWP3-W5E5 global meteorological data (available at https://data.isimip.org/search/tree/ISIMIP3a/InputData/climate/atmosphere/gswp3-w5e5/, last access: 9 May 2022) from 1980 to 2015 were used in all simulations in this study. The spatial resolution of the GSWP3-W5E5 data was 0.5^∘. Eight daily meteorological variables (downward shortwave radiation, downward longwave radiation, specific humidity, rainfall, snowfall, air pressure, wind speed, and air temperature) were used to run H08.

2.6 Reference yield data

To calibrate and validate the simulated crop yield, several yield datasets with different spatial resolutions were collected. The country-level yield data from FAO (available at https://www.fao.org/faostat/en/#data, last access: 9 May 2022) and grid-level (0.5^∘) yield data from the Global Dataset of Historical Yield (GDHYv1.2+v1.3) (Iizumi and Sakai, 2020) (available at https://doi.org/10.1594/PANGAEA.909132, Iizumi, 2019) for the period of 1986 to 2015 were used to evaluate model performance. FAO statistical yield was reported as fresh matter, whereas the model-simulated yield denotes dry matter. For consistency in the comparisons, as reported by Fader et al. (2010) and Müller et al. (2017), the FAO statistical yield was converted to dry matter with a crop-specific factor (e.g., 0.88, 0.88, 0.87, and 0.91 for maize, wheat, rice, and soybean, respectively) in accordance with Wirsenius (2000). The global dataset of historical yield for major crops (GDHY) is a spatially explicit dataset that converts the FAO annual national statistical yield to grid-level yield based on gridded net primary production estimated from several satellite products (Iizumi, 2020). The FAO statistical yield and GDHY yield provide valuable information for the evaluation of crop model performances at country and grid levels, respectively (Müller et al., 2017; Iizumi, 2020).

2.7 Simulation settings and yield processing

Individual simulations for maize, wheat, rice, and soybean were run under both rainfed and irrigation conditions from 1986 to 2015 on a daily scale. Details on simulation settings are listed in Table 1. The simulation was performed under the assumption that the four crops were planted and harvested in a hypothetical cropland of each grid cell. Under rainfed conditions, the crop growth was subjected to water stress; under irrigation conditions, there was no effect of water stress on crop growth. The yield processing was as follows.

First, the gridded yield (Yld) was aggregated from the simulated yield as follows:

$$\begin{array}{}\text{(18)}& \mathrm{Yld}={\displaystyle \frac{{\mathrm{Yld}}_{\mathrm{rain}}\times {\mathrm{Area}}_{\mathrm{rain}}+{\mathrm{Yld}}_{\mathrm{irri}}\times {\mathrm{Area}}_{\mathrm{irri}}}{{\mathrm{Area}}_{\mathrm{rain}}+{\mathrm{Area}}_{\mathrm{irri}}}},\end{array}$$

where Yld_rain and Yld_irri are the simulated yield under rainfed and irrigation conditions, respectively. Area_rain and Area_irri are the rainfed and irrigated harvest area per crop in a grid cell, respectively. The rainfed and irrigated harvest areas per crop were obtained from the MIRCA2000 dataset (Portmann et al., 2010) (available at https://www.uni-frankfurt.de/45218031/Data_download_center_for_MIRCA2000, last access: 9 May 2022).

Then, the national yield was aggregated from the gridded yield and weighted according to the crop-specific total harvest area. Because reference yield data have limited quality for marginal and small areas (Müller et al., 2017), we considered grid cells with a harvest area > 10 ha (Jägermeyr et al., 2021).

Finally, to ensure that the simulated data and reference data received similar treatment, we used the detrended yield when comparing time series variations in simulated yield and reference yields (Müller et al., 2017). In accordance with the methods of previous studies (Müller et al., 2017; Iizumi et al., 2013, 2014a), the moving average method was used to remove the trends. Specifically, similar to Müller et al. (2017), the anomaly yield was calculated by subtracting the moving average of a 5-year window.

3 Results and discussion

3.1 Effects of CO₂ fertilization and vapor pressure deficit

When only considering the CO₂ fertilization effect (simulation C), there was a positive impact on crop yield compared to default simulations (simulation D) (Fig. 1). In addition, similar to previous studies (e.g., Deryng et al., 2016), the CO₂ fertilization effect is larger for C₃ crops (wheat, rice, and soybean) than for C₄ crops (maize). In contrast, when only considering the vapor pressure deficit effect (simulation V), there was a negative impact on crop yield in comparison with default simulations. When considering the effects of both CO₂ fertilization and vapor pressure deficit, there was a positive impact on crop yield for the majority of the top 20 largest producer countries, while a negative impact was found for some countries (e.g., India and Egypt for maize). These impacts were also reflected in crop water productivity (CWP, defined as the ratio of crop yield to evapotranspiration). The averaged change in CWP in the top 20 largest producer countries was 4.8 %, −2.3 %, and 2.5 % for maize under simulations C, V, and CV, compared to simulation D (Fig. S2). The corresponding values were 6.4 %, −1.1 %, and 5.3 %; 5.8 %, −3.4 %, and 2.3 %; and 7.1 %, −3.6 %, and 3.4 % for wheat, rice, and soybean, respectively.

Figure 2 Comparison of mean simulated yield and mean FAO yield for the top 20 largest producer countries from 1986 to 2015. Dashed green and yellow lines indicate ± 10 % and ± 20 % differences, respectively. SIM denotes simulated yield, and FAO denotes reported yield from FAO. (a) Maize, (b) wheat, (c) rice, and (d) soybean.

Figure 3 Time series of detrended maize yield anomalies from simulation CVC (red), simulation D (blue), and FAO (green) for the top 20 largest producer countries. Y, yield; R, correlation coefficient; RMSE, root mean square error.

3.2 Model calibration

As shown in Figs. 1 and 2, the calibrated model (simulation CVC) showed better agreement with the FAO statistics of the mean national yield for the top 20 largest producer countries per crop (explaining approximately 88 %, 86 %, 93 %, and 99 % of global maize, wheat, rice, and soybean production, respectively). First, the mean bias (difference between mean national yield of simulation and mean national yield of FAO) of the 20 largest producer countries was considerably reduced to 2 %, 2 %, −2 %, and −1 % for maize, wheat, rice, and soybean, respectively. Second, the corresponding coefficient of determination (R²) values of the mean national yield of simulation and the mean national yield of FAO increased from 0.01 to 0.93, 0.21 to 1.00, 0.06 to 0.99, and 0.14 to 0.96 for maize, wheat, rice, and soybean, respectively. Third, the corresponding root mean square error (RMSE) decreased from 7.1 to 1.8, 2.2 to 0.3, 2.7 to 0.4, and 2.3 to 0.4 t ha⁻¹ for maize, wheat, rice, and soybean, respectively. These results suggested that the calibrated simulation could reliably reproduce the long-term averaged historical yield for the four major crops at the national level.

Figure 4 Comparison of R and RMSE values of time series of detrended yield anomalies between this study (SIM denotes simulation CVC) and Jägermeyr et al. (2021) (J21). Yellow bar denotes ensemble mean of different crop models used in the work of Jägermeyr et al. (2021). Error bars indicate maximum and minimum values among different crop models.

To investigate the capacity to reproduce the temporal variability in crop yield, time series of detrended yield anomalies in simulation data and FAO data for the top 20 largest producer countries per crop are presented in Fig. 3 for maize and Figs. S3–S5 for wheat, rice, and soybean, respectively. With regard to the ability to capture interannual variation in FAO yield, the model showed better performances for maize, wheat, and soybean than for rice. For example, positive correlations were found in 18, 16, 11, and 16 of the top 20 largest producer countries, with mean correlation coefficient (R) values of 0.48, 0.40, 0.31, and 0.37 for maize, wheat, rice, and soybean, respectively. The calibrated model showed better performance (increased R and decreased RMSE) than the default model in the majority of the 20 countries. Note that the calibrated model showed a similar performance to that of the default model in some countries (e.g., in USA, France, Ukraine, and Canada for maize) because the default simulations were already comparable to yield reported by the FAO, meaning that the calibration resulted in limited improvement (see Figs. 1a and 2a).

The R and RMSE values of time series of detrended yield anomalies between simulated yield and FAO yield for the top five largest producer countries per crop are summarized in Fig. 4. These countries were selected to make the data comparable with the latest global crop model intercomparison study by Jägermeyr et al. (2021), which includes 11 crop models for the period 1980–2010 (Fig. S10 in Jägermeyr et al., 2021). Overall, the R and RMSE values of our simulations were within the range of current mainstream crop models reported by Jägermeyr et al. (2021). For maize, wheat, and soybean, the R and RMSE values of our simulation were comparable with the ensemble means of different crop models reported by Jägermeyr et al. (2021); for rice, our simulation showed higher R values (except in Bangladesh and China) and lower RMSE values. However, the metric scores of our calibrated model and the other crop models in the work of Jägermeyr et al. (2021) remained low (e.g., few countries had R values > 0.5). This finding suggested that current crop models continue to experience difficulty in fully capturing the interannual variation in the historical yield because crop models only reflect the interannual climate signals in the simulated yields (Jägermeyr et al., 2021). This also indirectly implied that climate variation might not be the main driver of interannual yield variation for the major producer countries.

Figure 5 Relationship between maize yield (red: simulation CVC; blue: FAO) and total precipitation in the growing season from 1986 to 2015 for the top 20 largest producer countries. Y, yield; P, precipitation; R, correlation coefficient.

Figure 6 Relationship between maize yield (red: simulation CVC; blue: FAO) and mean air temperature in the growing season from 1986 to 2015 for the top 20 largest producer countries. Y, yield; T, air temperature; R, correlation coefficient.

To validate further the above assumption, we investigated the impacts of climate variables (i.e., precipitation and air temperature) on interannual yield variation by analyzing the correlations of total precipitation and mean air temperature in the growing season with the annual yield per crop. Using maize as an example (Fig. 5), there were no statistically significant relationships (p > 0.05) between precipitation and FAO statistical yield for most of the top 20 largest producer countries (17/20). Significant positive correlations between precipitation and the FAO statistical yield (p < 0.05) were found in only three countries: Romania, Hungary, and Serbia. The crop yield estimation relies on water availability; therefore, the variation in yield simulation largely reflects variation in precipitation. Accordingly, we observed good simulation performance in those three countries (Fig. 2a) with a clear correlation between FAO yield and precipitation (Fig. 5). Also, there were no statistically significant relationships between air temperature and FAO statistical yield for most of the top 20 largest producer countries (13/20) (Fig. 6). Similarly, there were no statistically significant correlations between precipitation and air temperature and FAO statistical yield in most countries for wheat, rice, and soybean (see Figs. S6–S11).

3.3 Comparison with GDHY gridded yield

Spatially explicit yield data enabled us to more fully evaluate the spatial distribution of model simulations. We compared the spatial distribution between simulated crop yield (simulations D and CVC) and the GDHY yield dataset. Using maize as an example, an apparent overestimation was detected in many parts of the world (e.g., China, Argentina, Brazil, India, Indonesia, Thailand, Mexico, and most countries in Africa) in the default simulations (Fig. 7a). In contrast, the calibrated simulation (Fig. 7b) showed a spatial pattern similar to the GDHY yield data (Fig. 7c). For the yields of wheat, rice, and soybean, the spatial distribution after improvement also showed a pattern similar to the GDHY yield data (Figs. S12–S14).

Figure 7 Spatial distribution of the mean (1986–2015) simulated yield of maize. (a) Simulation D; (b) simulation CVC; and (c) GDHY yield. Units in the legend are tonnes per hectare (t ha⁻¹).

In accordance with the method of Müller et al. (2017), we conducted a grid-level time series analysis of the correlations of the detrended yield between simulated and GDHY data (Fig. 8) to identify further the differences in the two yield datasets. Using maize as an example (Fig. 8a), statistically significant correlations (p < 0.1) were observed in a wide range of regions (e.g., northeastern USA, southern Europe, northeastern China, southern Brazil, eastern Argentina, southern Africa, and eastern Australia) (Fig. 8a), corresponding to 31 % of the total grid cells. Notably, there were also substantial differences in a considerable number of locations without statistically significant correlations (p > 0.1) (e.g., southeastern USA, western and central Asia, Brazil, and central Africa), corresponding to 69 % of the total grid cells (Fig. 8a). Similar characteristics were found for wheat, rice, and soybean (Fig. 8b–d).

Figure 8 Time series correlation between simulated yield (simulation CVC) and GDHY yield after trend removal using a 5-year moving average. Gray areas indicate no statistically significant correlation between the two datasets (p > 0.1), and white areas indicate no yield data for that crop in at least one of the two datasets. Panels show the determination coefficients for (a) maize, (b) wheat, (c) rice, and (d) soybean, respectively.

Such similarities or discrepancies between two yield datasets have been observed previously (see Fig. 9 in Müller et al., 2017). For example, there were statistically significant correlations (p < 0.1) and no statistically significant correlations (p > 0.1) between two datasets developed by Iizumi et al. (2014b; an earlier version of GDHY used in this study) and Ray et al. (2012) in a wide range of regions. Such comparisons can help in identifying considerable disagreements in global estimates of the spatial distribution of crop yield (Kim et al., 2021). Because it is difficult to determine whether one of these estimates is better than the others, the disagreement between our simulation and the GDHY data does not necessarily indicate that our simulation quality is low.

3.4 Limitations

Although crop yield simulations were improved, there were several limitations because of the assumptions, methods, and datasets used in this study. First, in accordance with the methods of previous studies (Müller et al., 2017; Jägermeyr et al., 2021), yield calculation and aggregation were conducted with the assumption that the irrigated harvest area and total harvest area per crop did not change throughout the study period; this assumption was based on data availability. However, these aspects do change over time. To overcome the problems associated with such an assumption, dynamic harvest area data at annual intervals, as generated by Mialyk et al. (2022), should be considered in future studies. Second, our calibration was conducted at the national scale in accordance with the method of Fader et al. (2010) rather than using finer spatial scales (e.g., subnational or grid-level), which increased the uncertainty of the yield simulations within each country. As shown in Fig. 7, the yield distribution is highly variable within a specific country. To incorporate the spatial heterogeneity in crop yield, parameter calibration should ideally be conducted at the grid-cell level (e.g., Iizumi et al., 2009; Sakurai et al., 2014). Although this approach has long-term promise, it is technically challenging because of uncertainty in the global gridded yield products and the potential for inflation in the parameter optimization calculation. In addition, the calibrated parameter reflected the mean average state, therefore potentially ignoring the yearly variation. Third, the reference dataset from GDHY does not represent purely observation-based yield, and, therefore, it is subject to errors or uncertainty resulting from its own methodology (e.g., errors in gross primary production and crop stress response) (Müller et al., 2017). Finally, our crop model is a simple model that does not fully represent the processes influencing crop growth. For example, we did not explicitly simulate N and P processes, although these effects are now reflected in the calibrated parameters (Fader et al., 2010). Additionally, the waterlogging effect is underrepresented in most crop models, including our model (Jägermeyr et al., 2021). Such physical mechanisms should be addressed in the development of future models.

4 Case study to estimate the contribution of irrigation to global food production

To demonstrate that the enhanced model can be applied for food–water nexus studies, we compared the predictions of a well-recognized study by Sibert and Döll (2010), which estimated the contribution of irrigation to global food production, with our predictions. This required a global crop yield model capable of estimating crop yield and explicitly dealing with the effect of irrigation.

Irrigation plays a critical role in global food production. The literature usually indicates that approximately 40 % of global total food production is from irrigated land (Postel et al., 2001; Siebert et al., 2005; Abdullah, 2006; Khan et al., 2006; Wada et al., 2013; Perrone, 2020; Ringler et al., 2020; Borsato et al., 2020), but the rationale and country-specific variation have not been fully explained. To our knowledge, Postel (1992) reported one of the first estimates, whereby approximately 36 % of the global food production was from irrigated land based on statistical data. Then, Siebert and Döll (2010) reported that irrigation contributed to approximately 33 % of the total global production. Here, we revisited the irrigation contributions for global production of maize, wheat, rice, and soybean using our improved model. Irrigation contribution in percentage (I) in a country (c) is defined as follows: I,c = $\frac{{\mathrm{Yld}}_{\mathrm{irri},c}\cdot {\mathrm{Area}}_{\mathrm{irri},c}}{{\mathrm{Yld}}_{\mathrm{irri},c}\cdot {\mathrm{Area}}_{\mathrm{irri},c}+{\mathrm{Yld}}_{\mathrm{rain},c}\cdot {\mathrm{Area}}_{\mathrm{rain},c}}\cdot \mathrm{100}\mathit{\%}$, where Yld_irri,c and Yld_rain,c are the irrigated and rainfed yields for a country, respectively; Area_irri,c and Area_rain,c are the total irrigated and rainfed harvest areas for a country, respectively.

Our results showed that the global average production levels from irrigated cropland were approximately 27 %, 30 %, 61 %, and 16 % for maize, wheat, rice, and soybean, respectively. These estimates were close to the estimates of Siebert and Döll (2010): 26 %, 37 %, 77 %, and 8 %, respectively. The similarities between these two studies mainly arose because both studies used data from Portmann et al. (2010) for crop-specific harvested area, and both models were calibrated with FAO data.

5 Conclusions

In this study, we determined the effects of CO₂ fertilization and vapor pressure deficit on crop yield using the global hydrological model H08. Then, we calibrated the yields of four major staple crops: maize, wheat, rice, and soybean. The calibrated national yield estimates generally showed good consistency with FAO statistical national yields. The calibrated grid-level yield estimates showed similarities in terms of spatial patterns and the reproduction of interannual variation compared with GDHY yield over a wide area, although there were substantial differences in other places. As reported in previous studies, the full reproduction of historical interannual yield variation remains challenging for global gridded crop modeling. Finally, we quantified the contributions of irrigation to the global production of maize, wheat, rice, and soybean; we explored the variations in irrigation contributions among countries. Our improvements provide a tool that can simultaneously simulate the water cycle and crop production while specifying irrigation water withdrawal in terms of the most detailed sources within a single framework, which will be beneficial for advancing global food–water nexus studies in the future (e.g., planetary boundaries, virtual water trade, and sustainable development goals).