Articles | Volume 13, issue 10
Model description paper
02 Oct 2020
Model description paper |  | 02 Oct 2020

MIROC-INTEG-LAND version 1: a global biogeochemical land surface model with human water management, crop growth, and land-use change

Tokuta Yokohata, Tsuguki Kinoshita, Gen Sakurai, Yadu Pokhrel, Akihiko Ito, Masashi Okada, Yusuke Satoh, Etsushi Kato, Tomoko Nitta, Shinichiro Fujimori, Farshid Felfelani, Yoshimitsu Masaki, Toshichika Iizumi, Motoki Nishimori, Naota Hanasaki, Kiyoshi Takahashi, Yoshiki Yamagata, and Seita Emori

Future changes in the climate system could have significant impacts on the natural environment and human activities, which in turn affect changes in the climate system. In the interaction between natural and human systems under climate change conditions, land use is one of the elements that play an essential role. On the one hand, future climate change will affect the availability of water and food, which may impact land-use change. On the other hand, human-induced land-use change can affect the climate system through biogeophysical and biogeochemical effects. To investigate these interrelationships, we developed MIROC-INTEG-LAND (MIROC INTEGrated LAND surface model version 1), an integrated model that combines the land surface component of global climate model MIROC (Model for Interdisciplinary Research on Climate) with water resources, crop production, land ecosystem, and land-use models. The most significant feature of MIROC-INTEG-LAND is that the land surface model that describes the processes of the energy and water balance, human water management, and crop growth incorporates a land use decision-making model based on economic activities. In MIROC-INTEG-LAND, spatially detailed information regarding water resources and crop yields is reflected in the prediction of future land-use change, which cannot be considered in the conventional integrated assessment models. In this paper, we introduce the details and interconnections of the submodels of MIROC-INTEG-LAND, compare historical simulations with observations, and identify various interactions between the submodels. By evaluating the historical simulation, we have confirmed that the model reproduces the observed states well. The future simulations indicate that changes in climate have significant impacts on crop yields, land use, and irrigation water demand. The newly developed MIROC-INTEG-LAND could be combined with atmospheric and ocean models to develop an integrated earth system model to simulate the interactions among coupled natural–human earth system components.

1 Introduction

The problems associated with climate change are related to the various processes involved in natural and human systems, as well as their interconnections. Changes in the climate system are caused by greenhouse gas emissions and changes in land use resulting from human activity (Collins et al., 2013). At the same time, climate change impacts natural and human systems in a variety of ways (e.g., Arent et al., 2014; Porter et al., 2014; Romero-Lankao et al., 2014). According to research on the linkage of various risks caused by climate change (e.g., Yokohata et al., 2019), changes in the climate system affect the natural environment, leading to changes in the socioeconomic system and finally impacting human lives.

One of the factors that plays an essential role in the interaction between the natural and human systems is land use (van Vuuren et al., 2012; Rounsevell et al., 2014; Lawrence et al., 2016). In general, changes in land use are driven by changes in various socioeconomic factors, such as an increase in food demand (Foley et al., 2011; Weinzettel et al., 2013; Alexander et al., 2015). At the same time, changes in the climate system affect the water resources available to agriculture and the size of the food supply through changes in crop yield (Rosenzweig et al., 2014; Liu et al., 2016; Pugh et al., 2016), significantly affecting human land use (Parry et al., 2004; Howden et al., 2007). Furthermore, climate mitigation measures often include the use of biofuel crops, which can significantly influence human land use (Smith et al., 2013; Humpenöder et al., 2015; Popp et al., 2017). On the other hand, land-use change is known to have biogeophysical and biogeochemical effects on the earth system (Mahmood et al., 2014; Chen and Dirmeyer, 2016; Smith et al., 2016), as changes in land use bring about changes in surface heat and water budgets, which, in turn, affect air temperature and precipitation (Feddema et al., 2005; Findell et al., 2017; Hirsch et al., 2018). Changes in land use also affect the terrestrial carbon budget, thereby influencing the concentration of greenhouse gases (GHGs) in the atmosphere (Brovkin et al., 2013; Lawrence et al., 2016; Le Quéré et al., 2018). It seems clear, then, that climate change induces land-use change by affecting various human activities, and that human land-use change affects changes in the climate system (Hibbard et al., 2010; van Vuuren et al., 2012; Alexander et al., 2017; Calvin and Bond-Lamberty 2018, Robinson et al., 2018).

Various numerical models have been developed to describe the interaction between natural and human systems in order to project future conditions as they relate to climate change (van Vuuren et al., 2012; Calvin and Bond-Lamberty 2018). Generally, in models dealing with the details of natural systems, elements related to human activity are simplified, and in models dealing with the details of human activities, elements related to natural systems tend to be likewise simplified (Müller-Hansen et al., 2017; Robinson et al., 2018). An earth system model (ESM) describes in detail the physical and carbon cycle processes in a natural system. A number of ESMs take human activities into consideration (Calvin and Bond-Lamberty 2018). The iESM project (Collins et al., 2015) is based on a CESM (Community Earth System Model Project, 2019) that incorporates GCAM (Calvin, 2011; Wise et al., 2014), an integrated assessment model (IAM) that provides a comprehensive description of human economic activities. With iESM, it is possible to capture the various interactions between the natural environment and human economic activities (Collins et al., 2015), but the model used to indicate the impact of climate change on water resources and crops is rather simplified (Thornton et al., 2017; Robinson et al., 2018; Calvin and Bond-Lamberty 2018).

IAMs consider supply and demand equations across the entire range of economic transactions and calculate the changes in surface air temperature resulting from increased GHGs in the atmosphere (Moss et al., 2010). IAMs can also project future changes in human land use (Wise and Calvin, 2011; Letourneau et al., 2012; Hasegawa et al., 2017). In general, however, IAMs simplify processes related to the natural environment (water resources, the ecosystem, crop growth, etc.) (Robinson et al., 2018) and thus do not explore the interactions between the natural and human systems on a spatially disaggregated basis (Alexander et al., 2018).

Many models for predicting changes in human land use have been developed (e.g., Hurtt et al., 2006; Lotze-Campen et al., 2008; Havlík et al., 2011; Wise and Calvin 2011; Meiyappan et al., 2014; Dietrich et al., 2019). Among these, the LPJ-GUESS and PLUMv2 coupled models are able to consider spatially specific interactions between changes in vegetation, irrigation, crop growth, and land use (Engström et al., 2016; Alexander et al., 2018). However, LPJ-GUESS (Olin et al., 2015) is a dynamic vegetation model that is incapable of exploring interactions related to physical processes, such as biogeophysical effects or future changes in water resources. On the other hand, LPJmL is a well-established global dynamical vegetation, hydrology, and crop growth model that can also consider the nitrogen and carbon cycles (Rolinski et al., 2018; von Bloh et al., 2018). The output of LPJmL (Bondeau et al., 2007), such as crop yield, land and/or water constraints, and vegetation and soil carbon, is used in the land-use model MAgPIE (Lotze-Campen et al., 2008; Popp et al., 2011; Dietrich et al., 2013; Kriegler and Lucht 2015; Dietrich et al., 2019). Although the gridded information of LPJmL is linked to MAgPIE (Alexander et al., 2018), the land-use change calculated by MAgPIE is not communicated to LPJmL (one-way coupling), making interactive calculations using the dynamic vegetation, hydrology, crop growth, and land-use models impossible.

In this study, we develop a global model that can evaluate the spatially detailed interactions between physical and biological processes, human water use, crop production, and land use related to economic activities. The model is based on the land surface component of global climate model MIROC (Model for Interdisciplinary Research on Climate, version 5.0; Watanabe et al., 2010), into which we have incorporated water resources, land ecosystem, crop growth, and land-use models. In the integrated model, which we call MIROC-INTEG-LAND (MIROC INTEGrated LAND surface model version 1), the budgets of energy, water, and carbon are determined by consistently considering the processes related to land surface physics, ecosystems, and human activities.

Section 2 in this paper explains the overall structure of MIROC-INTEG-LAND. The component models of MIROC-INTEG-LAND (climate, land ecosystem, water resource, crop growth, and land use), here called “submodels”, are described in detail in Sect. 3. Special attention is given to the land-use submodel, as it was specifically developed for inclusion into MIROC-INTEG-LAND and is expected to play a pivotal role. The other submodels – the climate, water resources, crop growth, and land ecosystem models – are based on models developed in the course of previous research. Section 3 outlines how the submodels used here differ from the original models. Section 4 explains the numerical procedure used to combine the submodels in the integrated model. Section 5 describes the data used for the various inputs and boundary conditions required to operate the integrated model. Section 6 verifies model reliability by comparing historical simulation results with various observational data. A summary of the results from simulations by MIROC-INTEG-LAND of future conditions and a discussion of the interactions between climate and water resources, crops, land use, and ecosystem are presented in Sect. 7. Finally, in Sect. 8, we discuss possible research themes regarding the interaction between natural and human systems that can be addressed using MIROC-INTEG-LAND.

Figure 1Relationships among variables in MIROC-INTEG-LAND. Components of the integrated model (submodels) are shown as colored boxes. Climate (land surface) and water resource components are represented with HiGWMAT (Pokhrel et al., 2012a), which is based on the land surface model MATSIRO (Nitta et al., 2014) in a global climate model MIROC (Watanabe et al., 2010). Land ecosystem and crop growth components are represented with VISIT (Ito and Inatomi 2012) and PRYSBI2 (Sakurai et al., 2014), respectively. The land-use model TeLMO (Terrestrial Land-use MOdel) is developed in this study. Inputs into the model are shown as boxes of climate and socioeconomic scenarios. Solid arrows between the boxes indicate the exchange of variables between the submodels. Dashed arrows indicate the input variables of the submodels.


2 Overall features of MIROC-INTEG-LAND

2.1 Model structure

The distinctive feature of MIROC-INTEG-LAND (Fig. 1) is that it couples human activity models to the land surface component of MIROC, a state-of-the-art global climate model (Watanabe et al., 2010). The MIROC series is a global atmosphere–land–ocean coupled global climate model, which is one of the models contributing to the Coupled Model Intercomparison Project (CMIP). MIROC's land surface component MATSIRO (Minimal Advanced Treatments of Surface Interaction and Runoff; Takata et al., 2003; Nitta et al., 2014) can consider the energy and water budgets consistently on the land grid with a spatial resolution of 1. MIROC-INTEG-LAND performs its calculations over the global land area only, and neither the atmosphere nor ocean components of MIROC are coupled. One of the advantages of running only the land surface model is that it can be used to assess the impacts of land on climate change, taking into account the uncertainties of future atmospheric projections.

Human activity models are included in MIROC-INTEG-LAND: HiGWMAT (Pokhrel et al., 2012b), which is a global land surface model with human water management modules, and PRYSBI2 (Sakurai et al., 2014), which is a global crop model. In HiGWMAT, models of human water regulation such as water withdrawals from rivers, dam operations, and irrigation (Hanasaki et al., 2006, 2008a, b; Pokhrel et al., 2012a, b) are incorporated into MATSIRO, the abovementioned global land surface model. In PRYSBI2, the growth and yield of four crops (wheat, maize, soybean, rice) are calculated. In addition, TeLMO (Terrestrial Land-use MOdel), a global land-use model developed for the present study, calculates the grid ratio of cropland (food and bioenergy crops), pasture, and forest (managed and unmanaged) as well as their transitions. The land-use transition matrix calculated by TeLMO is used in the process-based terrestrial ecosystem model VISIT (Vegetation Integrative SImulator for Trace gases; Ito and Inatomi 2012).

In MIROC-INTEG-LAND, various socioeconomic variables are given as the input data for future projections. For example, domestic and industrial water demand is used in HiGWMAT. The crop growth model PRYSBI2 uses future GDP projections in order to estimate the “technological factor” that represents crop yield increase due to technological improvement. The land-use model TeLMO uses future demand for food, bioenergy, pasture, and roundwood, as well as future GDP and population estimates. For future socioeconomic projections, we use the scenarios associated with shared socioeconomic pathways (SSPs; O'Neil et al., 2017) and representative concentration pathways (RCPs; van Vuuren et al., 2011). These are generated by an integrated assessment model: AIM/CGE (Asia-Pacific Integrated Model/Computable General Equilibrium; Fujimori et al., 2012, 2017b).

Interactions of the natural environment and human activities are evaluated through the exchange of variables in MIROC-INTEG-LAND (Fig. 1). The calculations in HiGWMAT are based on atmospheric variables (e.g., surface air temperature, humidity, wind, and precipitation) that serve as boundary conditions. The HiGWMAT model calculates the land surface and underground physical variables for three tiles (natural vegetation, rainfed cropland, and irrigated cropland) in each grid; a grid average is calculated by multiplying the areal weight of the three tiles. In HiGWMAT, water is taken from rivers or groundwater based on water demand (domestic, industrial, and agricultural). Agricultural demand is calculated endogenously in HiGWMAT, and withdrawn water is supplied to the irrigated cropland area, which modifies the soil moisture. The operation of dams and storage reservoirs also modifies the flow of the river. Using the soil moisture and temperature calculated in HiGWMAT, the crop model PRYSBI2 simulates crop growth and yield. PRYSBI2 also uses the same atmospheric variables that are used as input data in HiGWMAT.

The land-use model TeLMO uses the yield calculated by PRYSBI2. In TeLMO, the ratios of food plus bioenergy crop, pasture, and forest in each grid are calculated based on socioeconomic input variables such as the demand for food, bioenergy, pasture, and roundwood, as well as crop yield and ground slope. TeLMO also calculates the transition matrix of land usage (e.g., forest to cropland, cropland to pasture), which is passed to the terrestrial ecosystem model VISIT to evaluate the carbon cycle. The land uses calculated by TeLMO are also used as the grid ratios of natural vegetation and cropland area (rainfed and irrigated) in HiGWMAT.

2.2 Novelty of MIROC-INTEG-LAND

An important feature of MIROC-INTEG-LAND is that the land allocation model is coupled to the state-of-the-art land surface model, and that the impact of future climate and socioeconomic changes on water resources and land use can be considered consistently. In general, future land-use changes are often assessed by using an IAM. However, as mentioned earlier, IAMs are not grid based, but rather they divide the world into dozens of regions and describe the entirety of economic activity in these regions. Therefore, IAMs have a simplified description of the processes related to water resources and crop growth. In contrast, MIROC-INTEG-LAND provides capabilities to calculate complex physical processes over the land and considers the changes in water resources, taking into account human activities such as irrigation and reservoir operation. Furthermore, process-based crop models allow for an explicit and detailed consideration of growth processes of five different crops.

For the projection of future land use, IAMs usually (1) calculate the area of agricultural land by using yield information averaged over these regions based on the balance between supply and demand and (2) allocate the agricultural land by using a downscaling approach (e.g., Hasegawa et al., 2017). As pointed out in previous studies (Alexander et al., 2017), the problem with this method is that it does not allow for an explicit consideration of spatiotemporal information such as yield and production cost when determining land-use change. The Food Cropland Model in TeLMO addresses this issue by making it possible to consistently consider the spatiotemporal information such as crop yields and the balance between supply and demand when allocating the agricultural land, by using the Food Cropland Down-scale Module and the International Trade Module as explained in Appendix B.

As for the projection of future land-use change, TeLMO enables the calculation of future land-use change as an offline simulation, by using the crop yield data calculated in advance. On the other hand, crop yield depends on the water resource availability that is affected by the changes in soil physical processes due to future climate change, as well as the changes in irrigated cropland area caused by the increases in future food demands. MIROC-INTEG-LAND couples the models of physical land-surface processes, human water management, and crop growth processes with the land-use allocation model to consider these various interactions, as explained above.

3 Submodels

3.1 Global land surface model with human water management HiGWMAT

The HiGWMAT model (Pokhrel et al., 2015) is a global land surface model (LSM) that simulates surface and subsurface hydrologic processes considering both the natural and anthropogenic flow of water globally (1 in latitude and longitude). It incorporates human water management schemes (Pokhrel et al., 2012a, b) into the global LSM MATSIRO (Minimal Advanced Treatments of Surface Interaction and Runoff) (Takata et al., 2003). In MIROC-INTEG-LAND, HiGWMAT calculates the physical states (based on the changes in the energy and water budgets), including human water use and management. In HiGWMAT, the biophysical fluxes are updated after water use and management processes are simulated (Pokhrel et al., 2012a). Since our previous publications provide a detailed description of the MATSIRO model (Takata et al., 2003), groundwater scheme (Koirala et al., 2014), and the human impact representations (Pokhrel et al., 2012a, b, 2015), we include here only a brief overview of these models or schemes.

3.1.1 MATSIRO land surface model

MATSIRO (Takata et al., 2003; Nitta et al., 2014) was developed at The University of Tokyo and the National Institute for Environmental Studies in Japan as the land surface component of the MIROC (K-1 Model Developers, 2004; Watanabe et al., 2010) general circulation model (GCM) framework. MATSIRO estimates the exchange of energy, water vapor, and momentum between the land surface and the atmosphere on a physical basis. The effects of vegetation on the surface energy balance are calculated based on the multilayer canopy model of Watanabe (1994) and the photosynthesis-stomatal conductance model of Collatz et al. (1991) following the scheme in the SiB2 model (Sellers et al., 1996). The vertical movement of soil moisture is estimated by numerically solving the Richards equation (Richards, 1931) for soil layers in the unsaturated zone. The original version of MATSIRO (Takata et al., 2003) did not include an explicit representation of water table dynamics. To represent surface and subsurface runoff processes, a simplified TOPMODEL (Beven and Kirkby 1979; Stieglitz et al., 1997) is used. The surface heat balances are solved by an implicit scheme at the ground and canopy surfaces in the snow-free and snow-covered portions (i.e., four different surfaces within a grid cell) to determine ground surface and canopy temperature. The temperature of snow is prognosticated by using a thermal conduction equation, and the snow water equivalent (SWE) is prognosticated by using the mass balance equation considering snowfall, snowmelt, and freeze. The number of snow layers in each grid cell is determined from SWE. The albedo of snow in the model is varied using an aging factor (Wiscombe and Warren 1980) and in accordance with the time since the last snowfall and snow temperature, considering the densification, metamorphism, and soilage of the snow.

3.1.2 Human water management schemes

The original MATSIRO was enhanced by Pokhrel et al. (2012a, b) through the incorporation of a river-routing model and human water management schemes (i.e., irrigation, reservoir operation, water withdrawal, and environmental flow requirement). The irrigation scheme is based on the soil moisture deficit in the top 1 m (i.e., the root zone) of the soil column; that is, irrigation demand is estimated as the difference between the target soil moisture set for each crop type and the actual simulated soil moisture (Pokhrel et al., 2012b). Irrigation water is added as sprinkler irrigation on top of vegetation; a part of this is lost as evapotranspiration, and the rest returns back to the soil column. Subgrid variability of vegetation is represented by partitioning each grid cell into three tiles: natural vegetation, rainfed, and irrigated cropland.

The crop growth module for irrigation water is based on the H08 model (Hanasaki et al., 2008a, b), where the crop vegetation formulations and parameters are adopted from the Soil and Water Integrated Model (SWIM) (Krysanova et al., 1998). The crop growth module for irrigation water in HiGWMAT estimates the cropping period that is necessary to obtain mature and optimal total plant biomass for 18 different crop types. Irrigation is activated during the entire growing season but only for the irrigated portion of a grid cell using a tile approach (Pokhrel et al., 2012a).

Crop growth considered in the irrigation scheme is simulated within the HiGWMAT model using a crop growth module, which differs from the crop scheme in PRYSBI2 that simulates crop yields (Sect. 3.2). The reasons why different crop models are used to calculate irrigation water (HiGWMAT) and crop yields (PRYSBI2) are that (1) HiGWMAT has been used as a crop model based on SWIM (and it has been validated that the water withdrawal in various regions is consistent with the statistical data; Pokhrel et al., 2012b), and (2) PRYSBI2 has been used as a crop model based on SWAT, and crop yield in PRYSBI2 has been calibrated using the agricultural statistics; Sakurai et al. (2014). MIROC-INTEG-LAND uses different crop models to obtain realistic water withdrawal in HiGWMAT and to calculate realistic crop yields in PRYSBI2. The differences in the formulation between the crop models in PRYSBI2 and HiGWMAT are that the former uses more detailed crop modeling of the two-layer crop canopy, Farquhar photosynthetic CO2 assimilation, and the cropping period based on Sacks et al. (2010) (see details in Sect. A2), while the latter employs the simpler crop modeling of the single-layer crop canopy, radiation-use efficiency-type biomass accumulation, and the hypothetical planting date that gives the highest yield under the given weather conditions (Okada et al., 2015).

The reservoir operation and environmental flow requirement schemes are based on the H08 model (Hanasaki et al., 2008a, b). The reservoir operation scheme (Hanasaki et al., 2006) is integrated within the TRIP global river-routing model (Oki and Sud, 1998) to simulate reservoir storage and release for grid cells that contain reservoirs. The reservoir database is taken from Lehner et al. (2011). Large reservoirs having a storage capacity greater than 1 km3 are explicitly simulated; medium-sized reservoirs with a storage capacity ranging from 3×106 to 1×109 m3 (Hanasaki et al., 2010) are considered to be ponds holding water temporarily and releasing it entirely during the dry season. The withdrawal module extracts the total (domestic, industrial, and agricultural) water requirements: first from river channels and surface reservoirs and then from groundwater; the lower threshold of river discharge prescribed as the environmental flow requirement is considered when extracting water from river channels. While irrigation demand is simulated by the irrigation module, domestic and industrial water uses are prescribed based on the AQUASTAT database of the Food and Agricultural Organization (FAO; see Pokhrel et al., 2012b). We use the same prescribed values for domestic and industrial water uses in both historical and future simulations, as future projections of water withdrawal are not available.

3.2 Global crop growth model PRYSBI2

PRYSBI2 (Process-based Regional-scale crop Yield Simulator with Bayesian Inference 2) (version 2.2) is a semi-process-based global-scale crop growth model in which the daily biomass growth and resulting crop yield are calculated for the same grid cell as HiGWMAT (1 in latitude and longitude) (Sakurai et al., 2014). In MIROC-INTEG-LAND, PRYSIB2 is used to calculate crop yields. The target crops are maize, soybeans, wheat, and rice. Daily biomass growth is calculated using daily meteorological data (precipitation, temperature, wind speed, humidity, solar radiation, and atmospheric CO2 concentration) according to the photosynthetic rate calculated by a simple big leaf model (Monsi and Saeki, 1953) and the enzyme kinetics model developed by Farquhar et al. (1980). To determine the water stress, the soil moisture and temperature calculated by HiGWMAT (Sect. 3.1) are used. In PRYSBI2, the planting date is given by using the data of Sacks et al. (2010). The harvesting date is determined by when the crops accumulate their total number of heat units (THU) up to the threshold values. Crop yields for each year are calculated from the aboveground biomass and harvest indexes (Sect. A2).

The process of fertilizer input is not included in this model. Rather, parameters relating to technological factors that include the effect of fertilizer are set and input into the model (Sect. A7). We call this model a semi-process-based model, because some of the parameters, including the parameters relevant to technological factors, are statistically estimated using historical crop yield data (Iizumi et al., 2014) for each grid cell by the DREAM (DiffeRential Evolution Adaptive Metropolis) algorithm (Vrugt et al., 2009). The parameters were estimated by Markov chain Monte Carlo (MCMC) methods with 20 000 steps for each grid cell (Sakurai et al., 2014). The parameter values of the technological factors in future scenarios are estimated as a linear function of the gross domestic product (GDP) of each shared socioeconomic pathway (SSP) for each country (see details in Sect. A7).

In the original photosynthesis model by Farquhar et al. (1980), the photosynthesis rate is directly stimulated by the increase in CO2 concentration, which is called the CO2 fertilization effect. However, it is also known that the CO2 fertilization effect is downregulated by environmental limitations such as sink–source balance and nitrogen supply (Ainsworth and Long, 2005). In this model, the downregulation of the CO2 fertilization effect is described as a function of atmospheric CO2 concentration, in which the potential photosynthesis rate (maximum carboxylation rate of Rubisco and the potential rate of electron transport) gradually decreases according to the increase in CO2 concentration (see Sect. A6).

The crop model used in this study is an updated version (version 2.2) of the model described in Sakurai et al. (2014) (which gives a detailed description of PRYSBI2 version 2.0) and Müller et al. (2017) (which gives a brief description of version 2.1). The structure of the model is quite similar to versions 2.0 and 2.1. However, there are some parts of the version 2.2 structure that are slightly different. In Appendix A, we present a summary of the model and identify the elements that differ from the earlier versions.

3.3 Global land ecosystem model VISIT

The functions of the natural land ecosystem and their environmental responses are simulated by the submodel VISIT (Vegetation Integrative SImulator for Trace gases) (Ito, 2010; Ito et al., 2018). In MIROC-INTEG-LAND, VISIT is used to calculate the carbon and nitrogen cycles. VISIT is a process-based terrestrial biogeochemical model that simulates the atmosphere–land-surface exchange of greenhouse gases such as CO2 and CH4 and trace gases such as biogenic volatile organic compounds. Carbon, nitrogen, and associated water cycles are fully simulated in the model using ecophysiological relationships but in a simplified manner. The model operates at the global scale with a spatial resolution of 0.5×0.5. The ecosystem carbon cycle is simulated using a box-flow scheme composed of three plant carbon pools (leaf, stem, and root) and two soil carbon pools (litter and humus). Photosynthetic carbon acquisition is a function of the leaf area index, light absorptance, and photosynthetic capacity, which respond to temperature, ambient CO2, and humidity. Soil carbon dynamics are simplified by the litter and humus scheme but works well to simulate microbial decomposition and carbon storage. The model has two layers, i.e., natural vegetation and cropland, at each grid that are weighted by a land-cover fraction to obtain the total grid-based budget. Impacts of land-use change on the ecosystem carbon budget are taken into account using a simple scheme by McGuire et al. (2001) in which typical fractionation factors are applied to deforested biomass (e.g., immediate emission, 1-, 10-, and 100-year pools). The difference in carbon emissions from primary and secondary forests is included by using a different biomass density; regrowth of abandoned croplands is also simulated as the recovery of the mean biomass of the natural vegetation in the same grid. For brevity, croplands are categorized into three types (rice paddy, other C3 crops such as wheat, and C4 crops such as maize); the crop calendar and management practices such as fertilizer input are simulated within the VISIT model (i.e., independent of PRYSBI2) in a conventional manner. Planting and harvest dates are determined by monthly mean temperature; country-specific fertilizer inputs derived from the FAO country statistics (FAOSTAT; FAO, 2019) are used. In PRYSBI2, the effects of fertilizer are included in the technological factors, and crop yields are calibrated based on the technological factors, as described in Sects. 3.2 and A7. On the other hand, VISIT has been applied and validated at various scales from flux measurement sites to the global scale (e.g., Ito et al., 2017) based on the treatment of fertilizer input, as described above. The consistent treatment of fertilizer processes in PRYSBI2 and VISIT should be important future work.

3.4 Land-use model TeLMO

In the course of developing the integrated terrestrial model MIROC-INTEG-LAND, we developed the Terrestrial Land-use MOdel (TeLMO) for projecting global land use with a resolution of 0.5×0.5. TeLMO projects land use in each grid cell based on socioeconomic data such as demand for food and biofuel crops obtained from the AIM/CGE (Fujimori et al., 2012, 2017b). In MIROC-INTEG-LAND, TeLMO is used to estimate land-use change. For long-term projections, TeLMO assumes that there is a preferential order to land use by humans (i.e., urban, food cropland, bioenergy cropland, pasture land, and managed forests). That is, it assumes that land is used in the order of highest to lowest value added per unit area. After allocating land use in this manner, TeLMO calculates a transition matrix for each grid in order to evaluate the impact of land-use change on terrestrial ecosystems. Details of the five models comprising TeLMO – (1) the Food Cropland Model, (2) the Bioenergy Cropland Model, (3) the Pastureland Model, (4) the managed forest model, and (5) the land-use Transition Matrix Model – are explained in Appendix B.

Figure 2The numerical simulation procedure in MIROC-INTEG-LAND. The order of the numerical integration is (1) TeLMO, (2) HiGWMAT + PRYSIB2, and (3) VISIT as described in Sect. 4. Boxes indicate the submodels and data. For the submodels, the name and time step of the models are indicated in the boxes. In the “Data” box, the name of the variable saved as a file is indicated. In the “Input data” box, information regarding the input data is indicated.


4 Numerical procedure of model coupling

In MIROC-INTEG-LAND, submodels with different time steps are executed simultaneously by exchanging variables as shown in Fig. 1. The numerical procedure for exchanging variables between the submodels is shown in Fig. 2. Exchanging variables among submodels is accomplished in one of two ways: online coupling or offline coupling (Collins et al., 2015). In online coupling, the values calculated by a submodel are exchanged with other submodels via internal memory (i.e., the values calculated in one subroutine are passed directly to other subroutines). In offline coupling, the output of a particular submodel is written to a file; the other submodels then read the file as needed. The far-right “Data” box in Fig. 2 indicates the files used for saving submodel output data. The arrows show the exchanges that are made. The arrows between one submodel box and another indicate online coupling; those between a submodel box and the data box indicate offline coupling. The flow of submodel calculations is described below.

4.1 TeLMO

The land-use model TeLMO (Sect. 3.4) calculates the areal fraction of each land use within a grid (natural vegetation, cropland, pasture, etc.) and the transitions among them once a year, using the decadal average of crop yields calculated by PRYSBI2. The start year of TeLMO calculation is 2005. Since the exchange of variables is not so frequent, TeLMO is coupled to the other models via offline coupling (as shown in Fig. 2). That is, the output of TeLMO (grid fraction of land uses and transitions) is written to files, and the other submodels read the files as necessary. As shown in the figure, TeLMO reads the output files of PRYSBI2 (crop yields) for its calculations.


HiGWMAT (Sect. 3.1), the global land surface model that considers human water management, is used to calculate the physical states (surface and soil temperature and moisture, as well as energy and water fluxes) at hourly to daily time steps. The crop model PRYSBI2 (Sect. 3.2) is used to calculate crop yields at daily time steps using the soil moisture and temperature values generated by HiGWMAT. Since the exchange of variables between HiGWMAT and PRYSBI2 is very frequent (i.e., daily), these two submodels are joined through online coupling.

As shown in Fig. 2, in the future simulations, the MIROC-INTEG-LAND calculations start with TeLMO (TeLMO is switched off before 2004). After the output of TeLMO is written to files, the online-coupled HiGWMAT and PRYSBI2 make their calculations using the land-use grid ratio produced by TeLMO. Once the output of the HiGWMAT-PRYSBI2 combination is written to files, TeLMO again starts it calculations for the next year using the 10-year output. The exchange continues in this fashion.


As shown in Fig. 2, VISIT (Sect. 3.3), the terrestrial ecosystem model, calculates the carbon and nitrogen cycles using the output of the land-use model TeLMO. In MIROC-INTEG-LAND, no variable exchange between HiGWMAT-PRYSBI2 and VISIT is performed at this stage since the structures of these two submodels differ significantly. In the current version of MIROC-INTEG-LAND, we first calculate the TeLMO-HiGWMAT-PRYSBI2 calculations until the year 2100, and then perform the VISIT calculations from preindustrial time (including spin-up simulations) to the end of the 21st century by using the TeLMO output. (TeLMO is used only for the future period, and LUH (Land Use Harmonized) data are used for other periods.)

4.4 Model coupling

The proper choice of coupling method depends on the specific features of the variable exchange between submodels (Collins et al., 2015). One of the advantages of offline coupling is that the structure of the original model (e.g., the relationships between the main program and the subroutines) can be preserved, at least to some extent, in the coupling. This is not the case for online coupling. For example, for online coupling, either the main program of the original model needs to be modified in order for it to serve as a subroutine or a special program for connecting stand-alone models (i.e., a coupler) needs to be developed. In MIROC-INTEG, offline coupling is suitable for coupling TeLMO since the model structure of TeLMO is different from the other submodels (TeLMO solves equations with various spatial resolutions: global 30 s, 0.5, and 17 regions; see Appendix B for details) and data exchange occurs only once per year (so that the calculation cost for the input/output procedure can be minimized). On the other hand, online coupling is appropriate for connecting HiGWMAT and PRYSBI2, since the structure of the two submodels is similar (spatial resolution with a global 1 grid), and the exchange of variables is frequent (daily). In MIROC-INTEG, some of the subroutines of the original PRYSBI2 models that calculate the crop growth processes are called from HiGWMAT.

5 Experimental settings

Since MIROC-INTEG-LAND is based on a global land surface model, atmospheric boundary data (hereafter “forcing” data) are required to operate the model. The global land surface model with human water management HiGWMAT uses atmospheric temperature, humidity, wind, and surface precipitation as the forcing data to calculate the physical processes. In this study, we use forcing data from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) Fast Track (Hempel et al., 2013). In ISIMIP, historical and future climate simulations by five global climate models (GCMs) with bias correction are used as the distributed forcing data. The methodology of bias correction is described in Hempel et al. (2013). The five GCMs include GFDL-ES2M (Dunne et al., 2012), HadGEM2-ES (Jones et al., 2011), IPSL-CM5A-LR (Dufresne et al., 2012), Nor-ESM (Bentsen et al., 2013), and MIROC-ESM-CHEM (Watanabe et al., 2011). Uncertainties in the atmospheric predictions of the model can be considered by using the output data from the various GCMs. In ISIMIP data, correction for model bias is based on historical observations (Hempel et al., 2013). Thus, we can expect that over- and underestimation errors are removed (at least to some extent).

Since the time interval in the original ISIMIP data is daily and the time step in the land surface model HiGWMAT is subdaily, we generated 3-hourly data from the ISIMIP Fast Track daily data, based on the methods described in Debele et al. (2007) and Willet et al. (2007), where diurnal variations are generated based on the daily mean data.

In order to obtain a stable state of model variables, we performed spin-up simulations following the procedure defined in the ISIMIP Fast Track protocols. We first generated detrended 20-year data using 1951–1970 forcing data. The 20-year dataset was then replicated and assembled back-to-back to obtain an extended dataset. The order of years was reversed in every other copy of the 20-year block in order to minimize potential discontinuities in low-frequency variability. The time duration of the spin-up simulations was 400 years for the land surface model HiGWMAT and the crop growth model PRYSBI2 and 3000 years (repeated 100 times using the first 30 years detrended climate) for the terrestrial ecosystem model VISIT. The spin-up time of VISIT is longer than that for the other submodels, because it requires more time to reach a stable state, especially in the case of soil organic carbon.

After the spin-up simulations, we performed historical (1951–2005) and future (2006–2100) simulations based on the ISIMIP Fast Track protocols. For the future simulations, we used the forcing data of the five global climate models based on four RCPs (van Vuuren et al., 2011) – RCP2.6, 4.5, 6.0, and 8.5 – corresponding to radiative forcings of 2.6, 4.5, 6.0, and 8.5 W m−2 in the year 2100, respectively.

In the historical simulations of HiGWMAT, we used the land-use data (grid ratio of natural vegetation, rainfed, and irrigated cropland) provided by the Land Use Harmonized (LUH) project (LUHv2h; Lawrence et al., 2016); TeLMO was switched off. In the future simulations of HiGWMAT, the rainfed and irrigation cropland area is varied according to the output of TeLMO (Sect. 3.4). Since TeLMO projects the future total cropland area (irrigated plus rainfed), the future irrigated area is calculated by multiplying the grid irrigation ratio (irrigated  (rainfed + irrigated)) and the total cropland area calculated by TeLMO. The grid irrigation ratio is calculated by using the irrigated and rainfed cropland area determined by LUHv2h in 2005 and is fixed throughout the future simulation period. Although TeLMO also calculates the future bioenergy cropland area, we assume that bioenergy cropland is all rainfed.

TeLMO starts its calculations in 2005. As input data for TeLMO, we use the output variables based on the shared socioeconomic pathways (SSPs; O'Neil et al., 2017) calculated by an integrated assessment model (AIM/CGE; Fujimori et al., 2017b). In this study, we use outputs of the SSP2 scenario calculated by AIM/CGE (Fujimori et al., 2017b). Since the RCP8.5 scenario is not available in SSP2, we use the output of the baseline scenario by AIM/CGE for the calculation of RCP8.5. TeLMO uses future projections of GDP per capita, demand for food and bioenergy crops, pasture, and roundwood (Sect. 3.4, Appendix B). AIM/CGE calculates the aggregated transactions associated with the activities of economic actors; the energy system is represented in detail by dividing the globe into 17 regions (Fujimori et al., 2012).

The terrestrial ecosystem model VISIT is forced by the same ISIMIP forcing data used in HiGWMAT (Hempel et al., 2013). In the historical simulations, VISIT uses the historical land-use data from LUHv2h (Lawrence et al., 2016), as described above. In the VISIT future simulations, the output variables calculated by TeLMO, such as land use (cropland, pasture, forest) and the transition matrix describing transitions from one use to another (see Sect. 3.4 for details) are used as the forcing data.

It should be noted that the socioeconomic scenario that is used in climate forcing data by ISIMIP Fast Track (Hempel et al., 2013) does not match exactly the SSP scenarios (O'Neil et al., 2017), because the former is based on CMIP Phase 5 (CMIP5; Taylor et al., 2012) and the latter on CMIP Phase 6 (CMIP6; Eyring et al., 2016). This should not be a serious problem because the atmospheric processes are not coupled, and the radiative forcing (i.e., the RCP scenarios) used in ISIMIP Fast Track and the SSP scenarios is consistent. The ISIMIP phase 3 (ISIMIP3;, last access: 20 July 2020), which recently started distributing the climate forcing data, uses CMIP6 GCM simulations based on the SSP scenarios and is consistent with the present study.

Figure 3Comparison of historical terrestrial water storage (TWS) simulated by MIROC-INTEG-LAND with GRACE satellite data. For each river basin, the panel to the right shows the seasonal cycle. The GRACE data shown are the mean of the mass concentration products from two processing centers: CSR and JPL. Simulated results are the average of five climate model simulations. Grey shading indicates the uncertainty range shown by 1 standard deviation from the mean.


6 Historical simulations and comparisons with observations


Offline simulations from the original MATSIRO and HiGWMAT models have been extensively validated with ground- and satellite-based observations of various hydrologic fluxes and forms of storage (e.g., river discharge, irrigation water use, water table depth, and terrestrial water storage (TWS)) at varying spatial domains and temporal scales in numerous global-scale studies (Felfelani et al., 2017; Pokhrel et al., 2016, 2017, 2015, 2012a, b; Veldkamp et al., 2018; Zaherpour et al., 2018; Zhao et al., 2017). For completeness, we provide here a brief evaluation of TWS and irrigation simulations, since TWS is an indicator of overall water availability in a region and a primary determinant of terrestrial water fluxes (e.g., evapotranspiration (ET) and river discharge), and irrigation is an important component of the global freshwater systems that share the largest fraction of human water use globally (Hanasaki et al., 2008a; Pokhrel et al., 2016). Figure 3 plots the comparison of simulated TWS with observations by the Gravity Recovery and Climate Experiment (GRACE) satellite for the 2002–2005 period. The results shown are spatial averages over 18 major global river basins selected by considering a wide coverage of geographical and climate regions (Felfelani et al., 2017; Koirala et al., 2014). For the GRACE data, we use the mean of mass concentration (mascon) products from the Center for Space Research (CSR; Save et al., 2016) at the University of Texas at Austin and the Jet Propulsion Laboratory (JPL; Watkins et al., 2015; Wiese, 2016; Yuan et al., 2016) at the California Institute of Technology. It is evident from Fig. 3 that the model accurately captures the temporal variations as well as the seasonal cycle of TWS in most basins. Certain differences between model and GRACE can be seen in basins such as the Brahmaputra, Huang He, and Volga river basins, but such disagreements have been commonly reported in the literature owing to limitations in model parameterizations in simulating TWS components (e.g., the representation of snow physics and human activities) and inherent uncertainties in GRACE data (Felfelani et al., 2017; Scanlon et al., 2018; Chaudhari et al., 2019).

Figure 4Comparison of irrigation demands simulated by MIROC-INTEG-LAND (a) with the results from offline simulations using HiGWMAT and (b) forced by observed climate forcing data (Pokhrel et al., 2015) for 1×1 grids shown as the mean for the 1998–2002 period.

Figure 4 compares the irrigation water demand simulated by MIROC-INTEG-LAND with the results from offline HIGWMAT simulation obtained from Pokhrel et al. (2015), which is forced by the observed climate data. It is evident from this comparison that the broad spatial patterns seen in the offline simulations are clearly captured by MIROC-INTEG-LAND. Certain disagreements are, however, apparent. For example, MIROC-INTEG-LAND tends to overestimate irrigation demand over highly irrigated areas in the central United States, northwestern India, parts of Pakistan, and northern and eastern China, which is likely due to the drier and warmer climate simulated by MIROC (Watanabe et al., 2010) in these regions. The total global irrigation demand simulated by MIROC-INTEG-LAND is 1750 km3, which is greater than the 1238±67 km3 from the offline simulations but falls near the upper bound of estimates by various other global studies (see Table 1 in Pokhrel et al., 2015). The overestimation comes primarily from the highly irrigated regions noted above. Given that our meteorological forcing data are from GCM simulations, we consider our results for both TWS and irrigation demand to be acceptable.


Figure 5 shows historical simulation results for crop yield using ISIMIP forcing data as the baseline climate during the period from 1981 to 2005. The historical simulation results were compared with the gridded global dataset of historical yield (Iizumi et al., 2013; Iizumi, 2017), which is a hybrid of satellite-derived vegetation index data and FAOSTAT (FAO, 2019). The spatial aggregation to the country scale was conducted by using the harvested area (Monfreda et al., 2008). The area of wheat was separated into spring and winter wheat by using their production proportions (United States Department of Agriculture, 1994).

Figure 5Comparison of model estimation with reference data on average yield during the period 1981–2005 for the top 10 countries producing each crop. The box plots show the median and range of model results estimated from the five GCM outcomes. The main production countries were identified according to the country-based harvested area for each crop.


The results of the comparison of crop yields show that the simulated yields in most countries were underestimated to some degree (Fig. 5). Notably, using WATCH Forcing Data as the reference data in the bias correction for the ISIMIP dataset tends to underestimate solar radiation compared to the observation data (Iizumi et al., 2014; Famien et al., 2018), which in turn causes an underestimation of crop yields. The uncertainty of the projected yields as measured by the differences in outcomes for the five climate forcings was relatively small. The reason for this is that ISIMIP climate forcing data were bias corrected using the same historical weather dataset and the same method. For all crops, most of the correlations between the simulated and reported data were distributed along the 1 : 1 line. These results indicate that the model is capable of capturing the relative spatial difference of long-term average crop yield across countries.

Figure 6Comparison of latitudinal distribution of gross primary production (GPP) in 2000–2010 with upscaled flux measurements (Model-Tree Ensemble (MTE); Beer et al., 2010) and satellite observation (MODIS; Zhao et al., 2005).



The VISIT model captured the spatial and temporal patterns of terrestrial ecosystem productivity and carbon budget with satisfactory accuracy. Figure 6 shows the latitudinal distribution of gross primary production for the 2000–2010 period in comparison to upscaled flux measurements (Beer et al., 2010) and satellite observation (Zhao et al., 2005). High productivity in the humid tropics and low productivity in the arid middle latitudes and arid, cold high latitudes were effectively reproduced by the model simulation, although mean global total GPP was slightly higher than the observation (127.5 Pg C yr−1 by VISIT, 114.0 Pg C yr−1 by flux upscaling, and 121.7 Pg C yr−1 by satellite). Global carbon stocks in vegetation and soil organic matter were estimated as 499 and 1308 Pg C, respectively, in 2010; this is comparable to the contemporary synthesis (Ciais et al., 2013). Because of historical atmospheric CO2 rise, climate change, and land-use change, substantial changes in terrestrial ecosystem properties were simulated (not shown). As demonstrated by model validation and intercomparison studies, the VISIT model allows us to effectively capture the terrestrial ecosystem functions under changing environmental conditions.

6.4 TeLMO

In Fig. 7, the cropland area simulated by TeLMO in MIROC-INTEG-LAND is compared with the cropland area reported in FAOSTAT (FAO, 2019) and to the area simulated by AIM/CGE (Fujimori et al., 2017b), whose output of food demand and GDP per capita is used as input in TeLMO. Output of the TeLMO 0.5 grid data is aggregated by country to facilitate comparison with the FAOSTAT data. In order to also compare the TeLMO 0.5 grid data with the AIM/CGE cropland area, we used 0.5 downscaled land-use data based on the AIM/CGE calculation. (The methodology of downscaling is described in Fujimori et al., 2017a.) With the adjustment parameter Cj, the cropland area in TeLMO in 2005 is the same as that of LUH (Lawrence et al., 2016). As shown in Fig. 7, MIROC-INTEG-LAND roughly reproduces the cropland area by country shown in FAOSTAT (FAO, 2019). The differences in the five climate forcings given to MIROC-INTEG-LAND cause variance in crop yields, which in turn results in the variance in cropland area results shown in Fig. 7.

Figure 7Comparison of historical cropland area simulated by MIROC-INTEG-LAND (red), AIM/CGE (blue), and FAOSTAT (black), using the ratio of cropland area to total area. For MIROC-INTEG-LAND simulations, the cropland area results for the five different climate forcings are shown.


Figure 8Same as Fig. 7 but for the comparison of historical pasture area simulated by MIROC-INTEG-LAND (red), AIM/CGE (blue), and LUH (black), using the ratio of pasture area to total area.


In Russia, Brazil, and Australia, the recorded cropland area (i.e., FAOSTAT) is within the range of the MIROC-INTEG-LAND cropland area simulations using the different climate forcings. In Brazil and Russia, the variations in cropland area are mainly due to the difference in climate forcings. In the United States, the reported cropland area in FAOSTAT (FAO, 2019) is closely reproduced by MIROC-INTEG-LAND until around 2010; however, the declining trend of cropland area in the second half is not effectively reproduced. The reason for the overestimation seen here may be related to the underestimating of crop yield in PRYSBI2 (Sect. 6.3). The slight overestimation of the global cropland area trend (Fig. 7h) may stem from the same cause. Also, in China, although there is a declining trend of cropland area in MIROC-INTEG-LAND, in reality the cropland area remained nearly constant until 2014 and increased slightly thereafter. The increase of cropland area in China is considered to be influenced by policy, which is not considered in TeLMO.

In MIROC-INTEG-LAND, TeLMO uses the food demand and GDP per capita calculated by AIM/CGE under the socioeconomic scenario SSP2 (Fujimori et al., 2017b). Therefore, the difference between TeLMO and AIM/CGE is due to the difference in crop yield as well as the mechanism for the allocation of agricultural land. As explained in Sect. B1, TeLMO can consider the spatial distribution of crop yield when allocating agricultural land. On the other hand, in AIM/CGE, land-use change is calculated by aggregating crop yield information in the regions where the model calculation is performed (AIM/CGE divides the world into 17 regions). In large countries such as Australia, Brazil, and Russia, the allocation method in TeLMO shows good performance.

Figure 9Same as Fig. 7 but for the historical forest area simulated by MIROC-INTEG-LAND (red), AIM/CGE (blue), and FAO (black), using the ratio of forest area to total area.


Figure 10Time series of changes in the climate system based on the forcings of the five climate models. Results shown are for (a) surface air temperature (K); (b) soil moisture in the top 300 mm of the soil column (mm), shown as an anomaly from first 20-year average; (c) irrigation water supply (km3 yr−1). Thin curves indicate the global average of results for each of the five climate model forcings. Thick curves show the overall average of results based on the five forcings. The colors indicate RCP2.6 (blue), RCP4.5 (green), RCP6.0 (orange), and RCP8.5 (red).


Figure 8 shows a comparison of TeLMO, AIM, and LUH data for pasture. Unlike cropland, pastures are compared with LUH data, because there are no long-term global observation data. TeLMO calculates pasturelands such that the area matches that in the AIM for the AIM calculation domain (17 regions around the world). Because AIM treats China and the United States as one region, the results of TeLMO and AIM for China, the United States, and the globe are almost the same. On the other hand, in Australia, TeLMO is closer to LUH. Similarly, Fig. 9 shows a comparison between TeLMO, AIM, and FAO data of forest area. TeLMO refers to MODIS data and calculates forest area by taking into account deforestation and changes in crop area. Some differences can be seen between TeLMO and FAO, probably because TeLMO refers to MODIS and not to FAO; however, the differences are relatively small. Given that its performance is similar to that of AIM/CGE, the TeLMO submodel in MIROC-INTEG-LAND can be considered useful for future land-use prediction.

7 Future simulations and interaction of submodels

In the MIROC-INTEG-LAND future simulations, the RCP2.6, 4.5, 6.0, and 8.5 scenarios provided by ISIMIP1 (Hempel et al., 2013) serve as the climate scenario, while the output of AIM/CGE (demand for food and bioenergy crops, pasture, wood, etc.) according to the four RCPs under SSP2 (Fujimori et al., 2017b) serves as the socioeconomic scenario. The results in this section provide an understanding of the interactions between climate, water resources, crops, ecosystems, and land use that MIROC-INTEG-LAND accommodates.

Figure 11Time series of changes in crop yield (unit: t ha−1) based on the forcings of the five climate models. Results shown are for (a) winter wheat, (b) spring wheat, (c) maize, (d) soybean, (e) rice, and (f) grid maximum value for the five crop types. Thin curves indicate the global average of results for each of the five climate model forcings. Thick curves show the overall average of results based on the five forcings. The colors indicate RCP2.6 (blue), RCP4.5 (green), RCP6.0 (orange), and RCP8.5 (red).


Figure 10 shows the various time series related to climate system change. Figure 10a depicts the change in surface air temperature used as forcing data in MIROC-INTEG-LAND. It is displayed as the deviation from the average value of the 10-year period around the start year of the future simulations (2005). As shown in Fig. 10a, the increase in average global land surface air temperature in 2100 is approximately 6 C for RCP8.5, 3 C for RCP6.0, 2.5 C for RCP4.5, and 1 C for RCP2.6. Figure 10b shows the change in soil moisture calculated by MIROC-INTEG-LAND. Although the annual variation of soil moisture is considerable, the global land average soil moisture content tends to decrease in the 21st century. The reduction in soil moisture is largest in the RCP8.5 scenario, where the rise in surface air temperature is substantial. Results for the irrigation water supply are shown in Fig. 10c. As indicated in Sect. 3.1, water is supplied from rivers to the soil through irrigation until the ratio of soil moisture reaches a certain threshold. The irrigated area is calculated by multiplying the cropland area (as calculated by TeLMO) by the irrigation ratio, a fixed value corresponding to the ratio of irrigation cropland area to the total cropland area in 2005. Therefore, the changes in irrigation water supply in Fig. 10c reflect the changes in the irrigation area and the irrigation water supplied from rivers to the soil to compensate for the decrease in soil moisture. Although the global total cropland area increases in the first half of the 21st century (Fig. 12), in regions with a high irrigation ratio (e.g., India, China), cropland area decreases by the end of the century (Fig. 12). As a consequence, the irrigation area in MIROC-INTEG-LAND decreases, and, accordingly, the irrigation water supply also decreases, as shown in Fig. 10c.

Changes in crop yield calculated for the various future scenarios are shown in Fig. 11. The crop growth model PRYSBI2 in MIROC-INTEG-LAND can calculate the yields (t ha−1) of four crops (wheat, maize, soybean, rice), with a clear distinction between winter and spring wheat (meaning five crops in all). In Fig. 11f, the global average of the grid maximum yield value among the crops, which is used in the TeLMO calculation, is also shown. As described in Sect. 3.2, the future simulations by PRYSBI2 take into account the effects of climate change, as well as the CO2 fertilization effects due to rising greenhouse gas concentrations (Sect. A6) and the increase in technical coefficients due to future technological improvement (Sect. A7).

Figure 12Time series of changes in cropland area based on the forcings of the five climate models. The vertical axis is the cropland area as a fraction of total land area. The results are for (a) food cropland area and (b) food + bioenergy cropland area. Thin curves indicate the global average of results for each of the five climate model forcings. Thick curves show the overall average of results based on the five forcings. The colors indicate RCP2.6 (blue), RCP4.5 (green), RCP6.0 (orange), and RCP8.5 (red).


As shown in Fig. 11a–e, the yields of each of the crops rise over the first half of the 21st century. This is due to the CO2 fertilization effect and technological improvement. In general, the increase in yield is more significant in the high-GHG scenarios such as RCP8.5 than in the low-GHG scenarios such as RCP2.6. Such differences can be considered to be due to the fertilization effect and impact of climate change, since all the RCPs feature the same technological coefficient under the same SSP scenario (i.e., SSP2). On the other hand, in the latter half of the 21st century, the negative impact of climate change on crop yield is evident. In the RCP8.5 scenario, in particular, crop yields decline sharply. PRYSBI2 results show that the crop type most sensitive to climate change is maize: in 2100, the yield of maize under RCP2.6 is highest, while the yield of maize under RCP8.5 is lowest.

Figure 12a shows the change in the food cropland area calculated by TeLMO. As described in Sect. 3.4 and Appendix B, TeLMO uses the yield calculated by PRYSBI2 (grid maximum value as shown in Fig. 11f) and the food demand output of AIM/CGE. As shown in the Fig. 12a, crop area increases to meet the increase in food demand in the first half of the 21st century. Compared to other RCP scenarios during this time period, the RCP2.6 scenario requires more food cropland area, since the increase in crop yield is smaller in the RCP2.6 scenario. In the second half of the 21st century, the food cropland area tends to decrease as crop yield increases more than food demand. The decrease is smallest under RCP2.6 and largest under RCP6.0, and RCP8.5 actually requires an increase in food cropland area; as in this scenario, crop yields decline late in the century. Although there are differences among the results using the five different climate model forcings (the thin lines in Fig. 12a), using the average value lines (the thick lines in the figure) for comparison indicates that, by the end of the 21st century, the food cropland area is largest under RCP8.5.

Figure 12b shows the time series of the sum of food and bioenergy cropland area calculated by TeLMO. As described in Sect. 3.4, TeLMO calculates the distribution of the global bioenergy cropland area needed to meet the bioenergy demand calculated by AIM/CGE. It is known that the future bioenergy cropland area will change substantially depending on crop yield, and it should be noted that the setting in which crop yield is calculated can significantly affect the bioenergy cropland area (Kato and Yamagata, 2014). As shown in Fig. 12b, the bioenergy cropland area is significantly increased under RCP2.6 and RCP4.5. These climate scenarios require large areas of bioenergy crops for future climate mitigation. Although the food cropland area tends to decrease in the late 21th century (except in the RCP8.5 scenario), more cropland area will be needed if we consider both food cropland and bioenergy cropland.

Figure 13Spatial distribution of land-use change (units: a ratio of the grid box area). The results are for (a, b) food cropland area and (c, d) bioenergy cropland area. Average of the five climate projection-based simulations under (a, c) RCP2.6 and (b, d) RCP8.5 scenarios in the 2090s.

Figure 14Temporal change in global carbon stock in (a, b) vegetation biomass and (c, d) soil organic carbon, with (red) and without (green) land-use change under (a, c) RCP2.6 and (b, d) RCP8.5 scenarios. Thick lines show the median and light zones show the maximum to minimum range of the five climate projection-based simulations.


Figure 13 shows the global distribution of changes in food and bioenergy cropland areas, using the difference in 10-year averages around 2100 and 2005. As described in Fig. 12a, RCP2.6 tends to reduce the food cropland area in the latter half of the 21st century. Figure 13a and b show that the food cropland area decreases in Africa, India, and China. As is explained in Appendix B, TeLMO relies on the premise that the distribution of food cropland area is determined by changes in crop yield, food prices, wages (corresponding to changes in GDP per capita), and the demand for food. Thus the decreases in food cropland area shown in Fig. 13a and b are due to the increase in yield (meaning demand can be met with less cropland area) and the increase in GDP per capita (which means the population engaged in agriculture decreases due to development) in the SSP2 scenario. It should be noted that the change in cropland area at a particular grid is not determined solely by food production (the product of cropland area and crop yield) at that grid, as TeLMO considers the food trade among the 17 regions. As shown in Fig. 12 and noted earlier, the food cropland area will increase in the late 21st century in the RCP8.5 scenario. Accordingly, in comparison to the RCP2.6 scenario, the food cropland area in South America and central Africa increases in the RCP8.5 scenario.

Figure 15Spatial distribution of land-use-induced changes in terrestrial ecosystem carbon stock. Results are for (a, b) vegetation biomass and (c, d) soil carbon stock. Average of the five climate projection-based simulations under (a, c) RCP2.6 and (b, d) RCP8.5 scenarios in the 2090s.

As shown in Fig. 13, bioenergy cropland areas increase in various regions, especially in the RCP2.6 scenario. As discussed in Appendix B, TeLMO assumes that biofuel cropland is allocated based on the Agricultural Suitability Index (Eq. B14), which is a function of the yield and price of the bioenergy crop, GDP per capita, etc. At the same time, TeLMO also assumes that regions with high biodiversity are protected, and calculations are performed so as not to allocate biofuel cropland to the protected areas as shown in Fig. B2 (Wu et al., 2019). As a result, bioenergy cropland area is allocated to regions where the agricultural index is high – northwest and southern South America, central Africa, and Australia – but it cannot be allocated to protected areas such as the Amazon.

Figures 14 and 15 show the effects of changes in food and bioenergy cropland area on the terrestrial ecosystem calculated by VISIT in MIROC-INTEG-LAND. The impact of land-use change on terrestrial ecosystems is evaluated by comparing the calculation with and without considering the land-use change. The global time sequence (Fig. 14) shows that the changes in food and bioenergy cropland area have a significant impact on terrestrial ecosystems, especially in RCP2.6, where the aboveground biomass will decrease by approximately 50 Pg C (about 10 % of the present biomass stock) by 2100 due to deforestation for land-use conversion. The decrease in soil carbon after deforestation is much smaller than the decrease in aboveground biomass, as the carbon supply from crop residue compensates for the soil carbon loss. Consequently, this simulation implies that the impacts of land-use change occur heterogeneously and differ in their magnitude and direction between vegetation and soil. Figure 15 shows the global distribution of the effect of land-use change on aboveground biomass and soil carbon. The impact on aboveground biomass is projected to be greater in northwest South America, central Africa, northeast North America, and Australia, where the bioenergy cropland area is expanding. In these regions, even under the mitigation-oriented scenario, considerable declines in ecosystem structure and functions would occur, leading to deterioration, for example, of habitats for natural organisms, water holding capacity, and soil nutrients. Consequently, these functional degradations would degrade ecosystem services such as biodiversity, regulation, and provision. On the other hand, in Asia, the decrease in food cropland area tends to increase the aboveground biomass in both the RCP2.6 and RCP8.5 scenarios, possibly leading to the enhancement of aboveground biomass and thus ecosystem services.

Figure 16 shows the results of simulations to evaluate the effects of climate change on crop yield, land use, and water demand. In Fig. 16, the RCP8.5 simulations with climatic factors (temperature, water vapor, wind speed, soil moisture, soil temperature) and CO2 concentration fixed at 2006 (noCL + noFE), those with fixed climatic factors (noCL + FE), and those with variable climatic factors and CO2 concentrations (CL + FE) are compared. The CL + FE simulations are the same as the RCP8.5 results shown in Fig. 11f (crop yield), Fig. 12a (food cropland area), and Fig. 10c (irrigation demand).

As shown in Fig. 16a, the crop yield is significantly larger in the noCL + FE experiment than in the CL + FE experiment. This result indicates that climate change can significantly reduce crop yields. One of the reasons for the observed reduction in crop yield in the CL + FE experiment is that the growing season is shortened due to an increase in surface air temperature, which adversely affects crop growth (Sakurai et al., 2014). The impacts of climate change on crop growth increase with increasing temperature, and in 2100, crop yields in the CL + FE experiment are projected to decrease by approximately 60 % relative to the yields in the noCL + FE experiments.

As shown in Fig. 16a, the crop yield was much smaller in the noCL + noFE simulations than that in the CL + FE simulations. The reason for the yield in the noCL + noFE experiment being smaller than that in the CL + FE experiment is because the crop yield increases due to the CO2 fertilization effect in the latter. The increase in crop yield in the noCL + noFE experiment is due to technological developments (Sects. 3.2 and A7). Although there is a great deal of uncertainty regarding the treatment of CO2 fertilizer effects in crop models (Sakurai et al., 2014), the increase in crop yields due to the CO2 fertilizer effect is significant in the simulations of MIROC-INTEG-LAND.

Due to the changes in crop yields resulting from the changes in climate and fertilization effects, future cropland area and irrigation demand will also change significantly. In the CL + FE experiment, the food cropland area (Fig. 16b) and irrigation demands (Fig. 16c) become larger than those in the noCL + FE experiments because of the larger decrease in crop yields due to the impacts of climate change (Fig. 16a). On the other hand, the noCL + noFE experiment requires more food cropland area (Fig. 16b) and irrigation demand (Fig. 16c) compared to the CL + FE experiment because of the smaller increase in crop yields, mainly due to the absence of CO2 fertilization effects (Fig. 16a). In summary, the changes in climate and CO2 fertilization effects are expected to have marked impacts on crop yields, land use, and water demands in the future.

8 Implications and future research

With MIROC-INTEG-LAND, it is possible to calculate the interaction between climate, water resources, crops, land use, and ecosystems. The discussion in Sect. 7 suggests the type of feedback processes that can occur. As shown in Fig. 11, future climate change can affect crop yields. Especially under a scenario of large temperature increases (RCP8.5), crop yields will decrease in the latter half of the 21st century (Fig. 11). Here, the influence of the CO2 fertilization effect is also a very important factor affecting future changes in crop yields (Fig. 16a). Changes in crop yields due to climate change also have a large impact on cropland area (Figs. 12, 16b). Future cropland area may increase in response to an increase in food demand due to population growth, as well as due to increases in biofuel crop cultivation in response to global warming countermeasures. Such an increase in cropland area will cause a concomitant increase in water demand due to an increase in irrigated cropland area (Figs. 10, 16c). In addition, an increase in cropland area can affect carbon uptake in terrestrial ecosystems (Fig. 14). Increased human water use and changes in terrestrial carbon uptake can further affect the water, crop yields, and carbon budgets on the land surface. A real novelty of MIROC-INTEG-LAND is that the availability of both water and agricultural land can be consistently considered in conjunction with changes in climate conditions.

While this study showed only the results of the SSP2 scenario, in the SSP3 scenario, where the world is divided, the demand for food will be greater and more cropland area will be needed (O'Neill et al., 2017). Investigating the impacts of various natural and socioeconomic factors (climate, irrigation, fertilization effects, population, food demands, etc.) on land-use change and land ecosystems is an important future research direction as an extension of the present study.

Figure 16Time series of changes in (a) cropland yield (maximum across five crops in each grid, t ha−1), (b) food cropland area (a fraction of total land area), and irrigation demand (km3 yr−1) based on the forcings of the five climate models under the RCP8.5 scenario. Simulations with climatic factors and CO2 concentrations fixed at 2006 (light green, noCL + noFE), those with climatic factors fixed (cyan, noCL + FE), and those with varying climate and CO2 concentrations (red, CL + FE).


In addition to analyzing interactions, it is crucial to analyze the impacts of climate change and the effectiveness of countermeasures using MIROC-INTEG-LAND. The combined impacts of climate change on water resources, crops, land use, and ecosystems can be mitigated by enhancing various adaptation measures. For example, the use of water resources to control crop yield loss, changes in cropping calendars, and breeding can reduce the adverse effects of climate change on food and land use. With MIROC-INTEG-LAND, it is possible to assess the efficiency of adaptation measures designed to address the impacts of climate change on water resources, crops, land use, and ecosystems (Alexander et al., 2018). With consistent consideration of climate change, water resources, and land use, the competition between water, food, and bioenergy use can be analyzed (e.g., Smith et al., 2010). The model also provides useful insights into the trade-offs of biodiversity loss from land-use change and the benefits of climate mitigation.

MIROC-INTEG-LAND provides a way to integrate various human activity models based on the global climate model as shown in Sect. 4. This paper introduced illustrative simulation results produced by our application of MIROC-INTEG-LAND as a land surface model driven by meteorological forcing data. We plan to extend the model by enabling it to consider the physical processes and carbon and nitrogen cycles in the atmosphere and ocean. The MIROC community has developed MIROC-ES2L, an earth system model for CMIP6 (Hajima et al., 2020). By incorporating the water resource model (HiGWMAT), the crop growth model (PRYSBI2), and the land-use model (TeLMO) used in MIROC-INTEG-LAND into MIROC-ES2L, we are developing an integrated earth system model that we call MIROC-INTEG-ES. In MIROC-INTEG-ES, the interactions between the earth system and human activities are consistently considered. By using this integrated earth system model, the impact of land-use changes on the climate system, including biogeophysical and biogeochemical effects (Lawrence et al., 2016), can be more consistently investigated.

Appendix A: Description of crop model PRYSBI2 version 2.2

In the following description, we present a summary of the crop model used in MIROC-INTEG-LAND (PRYSBI2 version 2.2) and identify the elements that differ from the earlier versions (version 2.0, Sakurai et al., 2014, and version 2.1, Müller et al., 2017).

A1 Input data

As input data, PRYSIB2 version 2.2 uses the cropping period based on the planting and harvesting date by Sacks et al. (2010). Soil field capacity (Scholes and Brown de Colstoun, 2011) and atmospheric data (average, maximum, and minimum daily temperature; daily shortwave and longwave radiation; daily humidity; and CO2 concentration) are also used as input data. We use the same atmospheric data as HiGWMAT described in Sect. 5 (i.e., ISIMIP Fast Track data by Hempel et al., 2013).

A2 Growing period, maturity, and harvest

The time of seedling emergence after the planting date is determined by a parameter relevant to the average period between planting and emergence (lemerge). The period from emergence to maturity is determined by the total number of heat units (THU) (Neitsch et al., 2005). The crop is mature when THU is equal to a threshold value (thutotal), at which point it is harvested. THU thresholds were estimated for each grid by performing calibration between 1980 and 2006, so that harvest dates fit the data from Sacks et al. (2010). If future projections are performed using this threshold value, then the harvest date will deviate from Sacks et al. (2010) because of the temperature rise in future climates (i.e., harvest dates become earlier due to the increase in temperature). Using the biomass values obtained at the time of crop maturity, the yield is calculated as follows:

(A1) Yield = hi base BIO above ( maturity ) ,

where “Yield” is the crop yield (kg ha−1), hibase is the harvest index, and BIOabove(maturity) (kg ha−1) is the aboveground biomass at the time of crop maturity. Although the harvest index changes according to atmospheric CO2 concentration in version 2.0, in version 2.2, for simplicity, it is fixed.

A3 Photosynthesis

The photosynthesis processes in version 2.2 are the same as in the previous versions. The photosynthesis rate is calculated according to the daily meteorological data. The instantaneous global radiation and temperature at time (t) of the day are estimated from the daily global radiation and daily maximum and minimum temperature on a given day (td) according to the method described by Goudriaan and van Laar (1994). The amount of photosynthetically active radiation, PARt,td (MJ m−2 s−1), intercepted by the leaf at time t on a given day td, is calculated using Beer's law (Monsi and Saeki, 1953). We used the model described by Baldocchi (1994) to calculate the photosynthetic rate.

A4 Temperature stress

The equations for the effects of temperature on the maximum carboxylation rate of Rubisco and dark respiration rate are changed from those in version 2.0. The influence of temperature on the maximum carboxylation rate of Rubisco and the potential rate of electron transport is given as follows (Kaschuk et al., 2012; Medlyn et al., 2002):

(A2) C v c max ( t , t d ) = exp TM t , t d - 25 ep v c max 298 R ( 273 + TM t , t d ) ,

(A3) C j max ( t , t d ) = exp E j max TM t , t d - 25 298 R TM t , t d + 273 1 + exp 298 S j max - H j max 298 R 1 + exp TM t , t d + 273 S j max - H j max TM t , t d + 273 R ,

where Cvcmax(t,td) and Cjmax(t,td) represent the effect of temperature on the maximum carboxylation rate of Rubisco and the potential rate of electron transport, respectively; TMt,td is the air temperature (C) at time t on day td; epvcmax, Ejmax, Sjmax, and Hjmax are parameters that describe the shape of the curve (Kaschuk et al., 2012; Medlyn et al., 2002); and R is the universal gas constant (8.314 J mol−1 K−1).

The influence of temperature on the dark respiration of leaves is given as

(A4) C dark ( t , t d ) = exp TM t , t d - 25 ep rd 298 R ( 273 + TM t , t d ) ,

where Cdark(t,td) represents the effect of temperature on dark respiration at time t on day td, and eprd is the parameter that describes the shape of the curve (Kaschuk et al., 2012).

The maximum carboxylation rate of Rubisco, the potential rate of electron transport, and the dark respiration rate are modified by temperature effects:


where Vcmax(t,td) is the maximum carboxylation rate of Rubisco, Jmax(t,td) is the potential rate of electron transport, and vcmax and jmax are the potential maximum carboxylation rate and the potential rate of electron transport, respectively. Wstress(td) represents water stress, which is explained in Sect. A5. Θ is the compensation variable (0–1) that represents the discrepancy between the ideal photosynthetic potential and the actual one. ξV and ξJ are photosynthesis compensation variables that change according to CO2 concentration. These variables (Θ, ξV, and ξJ) are described in the following section. The dark respiration rate is calculated as follows:

(A7) R d ( t , t d ) = rd C dark t , t d v cmax ,

where Rd(t,td) is the dark respiration rate (µmol m−2 s−1), and rd is the leaf respiration factor (Collatz et al., 1991; Sellers et al., 1996a, b). The maintenance respiration and growth respiration are also considered. The formulations of the respiration models are also the same as those of the previous versions.

A5 Soil water balance and water stress

In PRYSBI2, the calculation of water stress follows the SWAT (Neitsch et al., 2005) algorithm. In SWAT, the daily water stress is calculated according to soil water, soil characteristics (field capacity and water content at saturation), root depth, and crop field evapotranspiration. PRYSBI2 uses the soil water calculated in HiGWMAT as explained in Sect. 3.2. The crop field evapotranspiration is calculated in SWAT according to the leaf area index.

A6 Correction of parameters according to CO2 concentration

The correction of parameters based on CO2 concentration is included in the model using the following equations:


where ξV and ξJ are photosynthesis compensation variables, drvcmax and drjmax describe the parameters, ca is atmospheric CO2 concentration (mol mol−1), and cbase is the baseline atmospheric CO2 concentration (mol mol−1). In this model, if drvcmax and drjmax>0, ξV and ξJ decrease linearly with increasing atmospheric CO2. If drvcmax and drjmax=0, ξV and ξJ do not depend on atmospheric CO2. In these equations, rmax1 and rmax2 are the respective asymptotic lines. rθ is the parameter that determines the curvature of the lines; we set rθ=0.99. The parameters drvcmax and drjmax are based on the results of Ainsworth and Long (2005).

A7 Time trend of the parameter relevant to agricultural management

When using historical yield data to calibrate model parameters, we need to consider temporal trends in the effects of non-climatic factors. Crop yield should improve from year to year because of agricultural factors, such as the decrease in harvest loss and the use of improved crop cultivars and pesticides. We, therefore, assumed the following linear trend in non-climatic effects when evaluating the long-term yield data:

(A14) Θ = θ base + θ trend Year - y base ,

where Θ is the compensation variable (0–1) that represents the discrepancy between the ideal photosynthetic potential and the actual one, which is used in Eqs. (A5) and (A6); θbase is the value of Θ in year ybase and must be calibrated for each cell of the grid; θtrend is the annual increase in Θ due to non-climatic factors (which also must be calibrated for each cell of the grid); “Year” is the year; and ybase is the criterion year (2006). In this study, we analyzed the relationship between θbase and GDP for each crop and used the estimated relationship for future prediction.

Appendix B: Description of land-use model TeLMO

B1 Food Cropland Model

For each grid, TeLMO first allocates the area for urban use; it then allocates the area for food cropland. For the allocation of the urban area, we use the Land Use Harmonization phase 2 future data that are used in Coupled Model Intercomparison Project Phase 6 (CMIP6) (LUH2f; Lawrence et al., 2016). It is generally expected that the food cropland area is determined by the balance between the supply and demand for food crops. The estimation of the supply potential of food crops requires the spatial distribution of crop production, which is related to the natural environment. On the other hand, the balance between the supply and demand for food crops is influenced by socioeconomic factors (e.g., populations, crop prices) related to international food trade. For this reason, TeLMO projects future land-use change by allowing the Food Cropland Down-scale Module (Sect. B1.1), which projects the global cropland distribution at a resolution of 0.5 by considering environmental factors, to interact with the International Trade Module (Sect. B1.2), which describes food supply and demand based on the General Equilibrium Model by dividing the world into 17 countries/regions. The primary objective of using TeLMO is to describe the long-term trend in land-use change and not the detailed year-to-year variations in land-use change. Therefore, we use 10-year average values as input to the model.

A major feature of TeLMO is that it does not project the local cropland distribution by unidirectionally downscaling the total cropland area for countries/regions obtained by integrated assessment models. This is because the total cropland area for each country/region depends on the local distribution of the cropland area. Therefore, TeLMO consistently treats the cropland distribution calculated by the Food Cropland Down-scale Module and the total cropland area for countries/regions obtained from the International Trade Module to project future land-use change. The Food Cropland Down-scale Module and International Trade Module are explained below.

B1.1 Food Cropland Down-scale Module

The Food Cropland Down-scale Module divides the earth into 0.5×0.5 (latitude × longitude) grid cells (hereinafter “0.5 cells”) and calculates the percentage of each cell occupied by cropland. The percentage of cropland is estimated by calculating the probability that each 30′′×30′′ grid cell (hereinafter “30′′ cell”) is used as cropland and averaging these probabilities over the entire 0.5 cell. A 30′′ cell allocated to urban use is not used for cropland. The probability ri of a given 30′′ cell being used as cropland is calculated as

(B1) r i = 1 1 + exp 1.228 + 0.237 ϕ i - 0.206 p k y j / w k C j ,

where ϕ is the slope, y is the yield per unit area (t ha−1), p is the price of food crops, w is the wage, and C is an adjustment parameter. The subscript i identifies the 30′′ cell, j identifies the 0.5 cell containing the ith grid cell, and k identifies the country/region containing the ith and jth grid cell. The definition of countries/regions is the same as that used in AIM/CGE (Fujimori et al., 2012, 2017b). Equation (B1) is formulated based on the fact that the cropland area is determined as a function of slope, crop price, yield, and the wages of farmers. The first term of Eq. (B1) is defined as the agricultural suitability index (ASI), which represents the relationship between cropland area and the explanatory variables. The adjustment parameter Cj is used to reproduce the cropland area of LUH (Lawrence et al., 2016) in the base year 2005 and to connect the future TeLMO projection with the historical simulation.

The ASI is derived from a logistic regression analysis using past statistical data. We use the global 0.5 MODIS cropland area (Friedl et al., 2010) as the objective variable, and the Global 30 Arc-Second Elevation (GTOPO30; Verdin and Greenlee 1996), the FAOSTAT food crop yield and price (FAO, 2019), and GDP per capita as the explanatory variables. GDP per capita rather than the wages of farmers is used for the reason indicated in the discussion of Eq. (B4) below. The logistic regression coefficient was derived from 23 000 data values that were randomly selected from the set of global 0.5 grids at year 2005. A comparison of the MODIS cropland areas and the calculated ASI values is shown in Fig. B1. The 23 000 randomly selected cropland area values were sorted in descending order and divided into 10 categories, and the average MODIS cropland area and the average ASI-based cropland area in each category were compared. As shown in Fig. B1, the values calculated by the logistic regression effectively reproduce the distribution of the MODIS cropland area data.

In the MIROC-INTEG simulations, GTOPO30 (Verdin and Greenlee, 1996) is used for the slope ϕi, and the food price pk and wage wk are obtained in the International Trade Module as explained in Sect. B1.2. PRYSBI2 results (1.0 resolution, Sect. 3.2), converted to a resolution of 0.5, are used for the yield yj. In TeLMO, total food cropland area is projected by using the maximum yield across the five cereal types (winter and spring wheat, maize, soybean, and rice). The reason for this formulation is explained in Sect. B1.2. yj in Eq. (B1) is calculated from the yields of the five cereal types by PRYSBI2. As discussed above, TeLMO is a model that evaluates the long-term trend in land-use change. Therefore, the crop yield and wage wk in Eq. (B1) is the average value of 10 years (using the data from the 1 year to the 10 years before the calculation year).

Figure B1Comparison of the global MODIS cropland area and the calculated area using the agricultural suitability index (ASI). Here, 23 000 randomly selected cropland area values are arranged in descending order and divided into 10 categories; the average value of MODIS (black) and ASI values calculated by TeLMO (red) in each category are compared. The horizontal axis is the higher percentile of cropland area data that is randomly selected from the global 0.5 grids at year 2005.


Figure B2Global distribution of areas protected from bioenergy production.

The 0.5 cell cropland area (Rj) is calculated by averaging the cropland probability in each of the 30′′ cells (ri) as follows:

(B2) R j = i J i r i J i ,

where Ji is the number of i cells (3600) in each 0.5 cell. The adjustment parameter Cj in Eq. (B1) is set so that the cropland area in the first year of calculation equals the data from LUH2f (Lawrence et al., 2016).

As explained above, the cropland distribution Rj projected at a spatial resolution of 0.5 by the Food Cropland Down-scale Module is used in calculations in the International Trade Module (Sect. B1.2).

B1.2 International Trade Module

Our model was developed by extending one of the simplest of the basic models – the Ricardian model. The Ricardian model is a one-production-factor (productivity per capita), two-country/two-commodity (food and non-food) model that attempts to describe the essence of free trade behavior based on the theory of comparative advantage. Because of its simple structure, the Ricardian model can be extended to a multi-country and multi-commodity model (Ejiri, 2008). In the International Trade Module, we extend the Ricardian model to be a multi-country (the entire world)/two-commodity (food and non-food) general equilibrium model. In addition, we account for decreasing returns in terms of production efficiency following the approach of Ejiri (2008). That is to say, we assume that agricultural production efficiency declines with increasing cropland area (and, conversely, that agricultural production efficiency increases as cropland area decreases). For this reason, industrial specialization, which has been pointed out as a problem of the Ricardian model, is unlikely to occur.

In order to construct a multi-country/two-commodity model, the subscript k was used to indicate country/region (the same 17 countries/regions defined in AIM/CGE), and subscripts 1 and 2 were added to indicate agricultural and nonagricultural sectors, respectively. The prices and wages in Eq. (B1) are those in the agricultural sector, which are represented by p1,k and w1,k, respectively.

First, wages in the agricultural sector, w1,k, are defined by using labor input and gross domestic production (GDP). In the International Trade Module, economic variables (e.g., food prices, wages, labor, and GDP) are described as the relative ratio to the base year (2005), which is the first year of calculation. Here, we assume that the total labor population ratio (relative to the base year) equals the total population ratio (relative to the base year).

(B3) l 1 , k + l 2 , k = L k ,

where l1,k, and l2,k are the labor input of the agricultural and nonagricultural sectors, respectively, and Lk is the total labor population (Murakami and Yamagata, 2019). GDP can then be described as total domestic income:

(B4) GDP k = w 1 , k l 1 + w 2 , k l 2 ,

where the value calculated by AIM/CGE is used for GDPk (units: USD). If we assume that the wage (ratio relative to the base year) for the nonagricultural sector is the same as that of the agricultural sector, the agricultural worker wage w1,k is calculated as

(B5) w 1 , k = GDP k l 1 , k + l 2 , k = GDP k L k .

In other words, it is assumed that the change in agricultural worker wage (relative to the base year) is equal to the change in per capita GDP. It is known that the employment rate have changed by a small percentage in the past. However, it is difficult to project the future changes in the employment rate, and thus the employment rate is assumed to be constant in the standard CGE (Computable General Equilibrium) models (e.g., Fujimori et al., 2012). Similarly, it is not easy to confirm the historical changes in wages for each country nor to estimate their future change; thus, similar to that for employment rate, the future changes in wages are usually kept constant in the CGE models (e.g., Fujimori et al., 2012). It should be noted that a small increase in employment rate (compared to the base year) can slightly decrease the wages as indicated in Eq. (B4), possibly leading to an increase in cropland area (Eq. B1).

Next, the price for agricultural sector p1,k is calculated using the multi-country/two-commodity general equilibrium model. The prices for agricultural and nonagricultural sectors are calculated using Eqs. (B5) and (B6), respectively:


where x1,k and x2,k are the production index in the agricultural and nonagricultural sectors, respectively. Here, the production index in the agricultural sector in region k (x1,k) can be calculated as the sum of the products of 0.5 crop yield yj and cropland area Rj using Eq. (7):

(B8) x 1 , k = j K j y j R j ,

where Kj indicates the number of 0.5 cells within the country/region k (3600). As described above, the cropland distribution Rj generated by the Food Cropland Down-scale Module (Sect. B1.1) is used in Eq. (B7). The domestic price p in Eqs. (B6) and (B7) is expressed in terms of the local currency unit (LCU). This is converted to the international price P (USD) using the exchange rate π (LCU per USD) in Eqs. (B8) and (B9):


The price p and production index x can then be connected using a relational equation for the trade budget as follows. Imposing the condition that the international budget for any country is zero results in Eq. (B10) for the international balance of payments:

(B11) p 1 , k x 1 , k - X 1 , k + p 2 , k x 2 , k - X 2 , k = 0 ,

where X1,k and X2,k are the demands for each good in each region. As described previously, the output generated by AIM/CGE based on the socioeconomic scenario is used for food demand X1,k. In this study, livestock feed demand is not included in X1,k. The international balance of payments as shown in Eq. (B10) consists of the current, capital, and financial accounts. The imbalance in the international budget corresponds to foreign exchange reserve. The foreign exchange reserve changes over periods longer than 10 years, but it is not possible to predict its future variation and thus is not considered in the standard CGE models (e.g., Ejiri, 2008). In the real world, if foreign exchange reserve increases, the amount of import goods tends to be decreased, because money is not used for them. Consequently, in food importing countries, food production tends to be increased, possibly leading to an increase in cropland area.

In addition, the price p and product index x can be related through Eq. (B11) by expressing economic growth in terms of GDP:

(B12) GDP k = P 1 , k x 1 , k + P 2 , k x 2 , k .

In Eq. (B3) and Eqs. (B5)–(B11) above, the eight unknown values are p1,k, p2,k, x1,k, x2,k, l1,k, l2,k, πk, and X2,k. Of these, because the reference for the international price P is the United States (region index k=1), P1,1 and P2,1 (along with p1,1, p2,1) cannot be set. For this reason, the condition is imposed that total global net exports and imports are equal to zero:


As explained above, TeLMO uses 10-year averages as input to the model to represent long-term trends in land-use change (Sect. B1.1). We assumed that the global total production is equal to consumption, i.e., the total global net exports and imports equal to zero. In reality, there are certainly stock changes in various goods, but it would not be counterfactual to assume that they are net zero at longer timescale. The unknown values for p1,k, p2,k, x1,k, x2,k, l1,k, l2,k, πk, and X2,k are calculated by simultaneously solving eight equations, Eq. (B3) and Eqs. (B5)–(B11), for all 17 regions (k=1-17) subject to the conditions imposed by Eqs. (B12) and (B13). The p1,k, and w1,k values obtained from Eq. (B4) are entered into Eq. (B1). Finally, the share of cropland for each 0.5 cell Rj can then be calculated using Eq. (B2).

As explained in Sect. B1.1, TeLMO uses the maximum yield of five cereal types to project the total cropland area. Alternatively, it is possible to increase the number of agricultural sectors in Eqs. (B3)–(B12), solve the prices for each crops, and allocate the cropland area according to the ASI values for each crop. Although we attempted this formulation in the course of our development of TeLMO, it was found that the results were similar to those obtained from the current formulation. On the other hand, the solution of general equilibrium models did not converge in some cases, because the number of sectors increases in the equations. For this reason, we decided to adopt the current formulation while recognizing that calculating cropland areas for each crop is an important future work.

B2 Bioenergy Cropland Model

The Bioenergy Cropland Model uses 30′′ cells that are not assigned to urban use or food cropland use. Whereas adjustment parameter Cj in the Food Cropland Model (Eq. B1) could be set using observed cropland area for the first year of the TeLMO calculation (the base year 2005), there is no corresponding adjustment parameter in the case of bioenergy cropland, because sufficient cropland devoted to biofuel crops did not exist in the base year. Accordingly, the Bioenergy Cropland Model allocates bioenergy cropland around the globe so that the global total biofuel crop production equals the global total biofuel crop demand obtained by AIM/CGE. The Bioenergy Cropland Model uses the same formularization to that in the Food Cropland Down-scale Module (Sect. B1.1) to evaluate the probability of bioenergy cropland in 30′′ cells using the following equation:

(B15) r bio , i = C bio 1 + exp 1.228 + 0.237 ϕ i - 0.206 p bio , k y bio , j / w 1 , k ,

where ϕi is the slope in 30′′ cell i, pbio,k is the biofuel crop price in region k, ybio,j is the yield (t ha−1) of biofuel crops in 0.5 cells, and w1,k is the agricultural sector wage in region k. For the biofuel crop price pbio,k, the values generated by AIM/CGE are used. For biofuel crop yield ybio,j, the yield for miscanthus or switchgrass, whichever is greater in a given cell, is calculated for the entire globe by using the biofuel crop model developed in Kato and Yamagata (2014). The biofuel crop model in Kato and Yamagata (2014) considers the future changes in climate based on the RCP scenarios. In this study, we also consider the future changes in fertilizer input based on the SSPs adopted in Mori et al. (2018). Because of the uncertainty in future fertilizer application for crop management, we set the high end of the N fertilizer input threshold according to Tilman et al. (2011). The nitrogen fertilizer application was set to increase from the current level according to the increasing rate of GDP in the SSP2 scenario up to 160 kg(N) ha−1 yr−1 if the fertilizer input at the country level was below 160 kg(N) ha−1 yr−1 in the 2000s. Also, the phosphorus fertilizer input in each country was set to follow the same annual increase rate as the nitrogen fertilizer application.

Our use of the same formularization for the Food Cropland Model and the Bioenergy Cropland Model is based on the assumption that the factors determining both cropland areas are similar.

The adjustment parameter Cbio is set so that the global total biofuel crop production volume (product of yield and cropland area) equals the global total biofuel crop demand calculated by AIM/CGE:

(B16) k K all X bio , k = j J all y bio , j R bio , j ,

where Xbio,k is the biofuel crop demand for region k calculated by AIM/CGE and Kall and Jall are the total number of regions (17) and the total number of 0.5 cells (259 200), respectively. Rbio,j is the average percentage of bioenergy cropland for all 30′′ cells in a given 0.5 cell, where the individual 30′′ cell percentages are determined by Eq. (B14).

If bioenergy cropland were allocated based on the principle described above, a massive development of bioenergy cropland would occur in regions with high ecosystem production such as the Amazon. For this reason, the model accounts for protected areas that cannot be allocated as bioenergy cropland as shown in Fig. B2. Two sources were used for protected areas (Wu et al., 2019): the World Database for Protected Areas (WDPA) (IUCN and UNEP-WCMC, 2018) and the World Database of Key Biodiversity Areas (KBA) (BirdLife International, 2017). As of 2018, the WDPA covered an area of 33.6×106 km2, and the KBA covered an area of 19.9×106 km2. In this study, we did not consider the protected area for the calculation of the food cropland and pasture under the assumption that food has a higher priority than ecosystem protection.

B3 Pastureland Model

Whereas the Food Cropland Model uses statistical relationships between cropland area, yield, and economic variables, because reliable statistical data do not exist for pastureland, a simpler approach is taken to estimate pastureland. The probability of pastureland in each 30′′ cell is determined based on net primary production (NPP) and slope, which is given by

(B17) r past , i = C past , j × NPP j 1 + ϕ / 20 .

The denominator in Eq. (B16) reflects the fact that the use of land as pasture decreases with the angle of inclination, as is shown in the LUH2f data (Lawrence et al., 2016). The results of an offline simulation by VISIT (Ito and Inatomi 2012), assuming the entire world to be grassland, are used here for NPPj. The boundary condition of the VISIT offline simulations is fixed at year 2005. Cpast,j is the adjustment parameter for 0.5 cells. The value of Cpast,j changes from year to year. The adjustment parameter for the base year, Cpast,j(t=0), is set so that the pastureland distribution equals that of LUH2f (Lawrence et al., 2016) for the base year (2005). Adjustment parameters for years other than the base year, Cpast,j(t), are set by applying a proportionality factor, α(t), to the base-year parameter:

(B18) C past , j ( t ) = α ( t ) × C past , j ( t = 0 ) ,

where α(t) is set so that regional total pastureland area equals the regional total pastureland demand calculated by AIM/CGE. In other words, α(t) is set so that the condition

(B19) S past , k ( t ) = j J k R past , j ( t )

is met, where Spast,k(t) is the pastureland demand calculated by AIM/CGE for region k; Rpast,j(t) is the average of percentage of pastureland for all 30′′ cells (from Eq. B16) in a given 0.5 cell; and Jk is the total number of 0.5 cells in each region k.

B4 Managed Forest Model

In the Managed Forest Model, satellite data are used to determine forest area; the share of forest area where timber harvesting occurs is allocated as managed forest in the manner described below. The distribution of managed forests in 0.5 cells, Rmfr,j(t), is formularized in terms of the area of managed forests in the base year and the population density:

(B20) R mfr , j = A fr , j × ρ j * C manfr , k + ρ j * ,

where Afr,j is the area of managed forest in the 0.5 cells in the base year (2005) and ρj* is the mean population density in the 5×5 grid (2.5 cell) of cells centered on the 0.5 cell in question. Larger 2.5 cells were used instead of 0.5 cells based on the assumption that harvested timber is transported within an approximately 100 km radius and that the amount of harvested timber is determined by the population density in each 2.5 cell. The 100 km radius is estimated from the distance where the transportation cost of timber ( USD 1 km−1 t−1) is balanced with the price of timber ( USD 100 t−1). Here, the transportation cost and price of timber are estimated using the FAOSTAT data (FAO, 2019). Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data (Friedl et al., 2010) are used for the base-year forest area (2005), and data from Murakami and Yamagata (2019) are used for the population density (ρj*). Cmfr,k is an adjustment parameter that is set for each of the 17 regions (k) so that the managed forest area conforms to the roundwood demand Xmfr,k (kg yr−1) calculated by AIM/CGE. We use the region-level adjustment factors for managed forest (Cmfr,k), because the grid-level reference data are not available. In other words, Cmfr,k is set so that the total regional amount of harvested timber equals the regional total roundwood demand:

(B21) X mfr , k = j J k R mfr , j × B j L j ,

where Bj is the distribution of forest biomass (kg m−2) in 0.5 cells, calculated by VISIT (Ito and Inatomi, 2012) offline simulations assuming the entire world to be forest with the fixed boundary conditions (2005). Jk is the total number of 0.5 cells in each region k. Lj is the harvesting period (years), which is estimated as follows, based on theNPPj for 0.5 cells obtained from VISIT (Ito and Inatomi, 2012).

(B22) L j = NPP j < 4 500 / NPP j 4 NPP j 25 20 25 < NPP j

Lj reflects the fact that the harvesting period decreases with increases in net primary production, as is shown in the LUH2v data (Lawrence et al., 2016). The amount of forest harvested in a given year can also be calculated as Rmfr,j×Bj/Lj (kg yr−1) based on the distribution of managed forests Rmfr,j, forest biomass Bj, and the felling period Lj for 0.5 cells.

B5 Formulation of the Transition Matrix Model

Evaluating the impact of land-use change on terrestrial ecosystems requires not only the spatial distribution of land use but also information on the land-use transition. For example, in areas where shifting cultivation is practiced, a particular area may be cleared as cropland while another area is abandoned even though the overall cropland area within a cell does not change. In such cases, there is a transition from cropland to secondary land, which impacts the aboveground biomass and carbon budget. Thus, matrix information regarding the transition from one land use to another land use is essential.

For the land-cover types used in the transition matrix, we use the five classes (urban, cropland, pasture, secondary/primary land) used in the VISIT terrestrial ecosystem model (Ito and Inatomi, 2012). TeLMO forecasts eight land-cover types, including the previously described urban, cropland (food and bioenergy), pasture, managed forest, and unmanaged forest classes as well as “grassland” (obtained from MODIS satellite data, Friedl et al., 2010) and “other” land-cover types that are not used by humans (e.g., glaciers, lakes, and marshes, as defined by MODIS satellite data; Friedl et al., 2010). The correspondence between the land-cover types used in TeLMO and those used in the land-use transition matrix is presented in Table B1.

Table B1Correspondence of land-cover type in land-use model and transition matrix.

Download Print Version | Download XLSX

The primary/secondary land classes in the land-use Transition Matrix Model are defined as land that has never been used by humans or land that has been used at least once by humans, respectively. Here, unmanaged forest and grassland are classified as primary or secondary land based on data from LUH2f supplied by LUH2v (Lawrence et al., 2016). Unmanaged forest or grassland areas that are classified as secondary land in the base year (2005) remain classified as secondary land in subsequent years. In the case in which unmanaged forest or grassland area is classified as primary land in the base year, it is classified as secondary land if the area is converted to cropland or pasture and then later returned to being unmanaged forest or grassland. In TeLMO, land classified as other is considered the land that cannot be used by humans and is therefore not included in the land-use transition matrices.

The method used to create the land-use transition matrices is shown in Fig. B3. As explained above, TeLMO assumes that land is used in order of highest to lowest value added per unit area (i.e., urban, food cropland, bioenergy cropland, pastureland, managed forest, and unmanaged forest). Aligning these land-use classes with corresponding classes in the transition matrix (Table B1), the preferential order of the latter becomes urban, cropland (food + bioenergy), pasture, secondary land, and primary land. To calculate land-use transition matrices, the percent areas of the different land-cover types in each 0.5 cell in a given year are first sorted in order of preference (“Pre” in Fig. B3). In Fig. B3, the length of each colored bar represents the percent area of a given land-cover type. The sum of the percent areas for all land-use classes is 100 %. Next, the percent areas of different land-cover types in each 0.5 cell in the following year are again sorted in order of preference (“Post” in Fig. B3).

Figure B3Schematic diagram of land-cover transition. Details are explained in the main text.


As shown in Fig. B3, the percent areas of transitioned land defined in transition matrices can be calculated by comparing the percent areas for each land-cover type in a given year and the next year. For example, the area indicated in column “a” in Fig. B3 corresponds to the percent area of land that transitioned from pasture to cropland. Similarly, the area indicated in column “b” in Fig. B3 corresponds to the percent area of land that transitioned from secondary land to pasture. In this manner, it is possible to calculate the transition between land-cover types by assuming a preferential order to land use.

Shifting cultivation is taken into account when making the land-use transition matrices. We assume that the share of cultivated land does not change over time on the larger (i.e., 0.5 cell) scale. Data from Butler (1980) are used for the global allocation of shifting cultivation on this larger scale. Furthermore, in regions where shifting cultivation is practiced, we assume that cropland is used sequentially with a fixed rotation (Butler, 1980). Under this assumption, in areas where shifting cultivation is practiced, 1∕15 of the cropland area is newly cultivated, and 1∕15 of the cropland area is abandoned each year. Thus, 1∕15 of the cropland area is transitioned from secondary land to cropland, and 1∕15 of the cropland area is transitioned from cropland to secondary land. These transitions are added to the transition matrices for areas where shifting cultivation is practiced.

Code and data availability

The MIROC-INTEG-LAND source code for this study is available to those who conduct collaborative research with the model users under license from the copyright holders. For further information on how to obtain the code, please contact the corresponding author. The data from the model simulations and observations used in the analyses are available from the corresponding author upon request.

Author contributions

TY was responsible for the development and description of MIROC-INTEG-LAND and TeLMO. TY and AI carried out the experimental study and analyzed the results of numerical experiments of historical and future simulations. TK was responsible for the development and description of TeLMO. GS, MN, and TI were responsible for the development and description of PRYSBI2. AI was responsible for the development and description of VISIT. MO analysed the output of PRYSIB2 historical simulations. YS calculated the forcing data. EK carried out the simulations for the bioenergy crop yields. YP, FF, TN, YM, and NH were responsible for the development and description of HiGWMAT and analyzed the output of HiGWMAT historical simulations. SF and KT were responsible for the input data by AIM/CGE. YY, SE, and TY proposed the development of MIROC-INTEG-LAND and led the project of model development. All authors have read and approved the final article.

Competing interests

The authors declare that they have no conflict of interest.


This research is supported by the “Integrated Research Program for Advancing Climate Models (TOUGOU Program)” sponsored by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan. It was carried out as part of the Integrated Climate Assessment–Risks, Uncertainties, and Society (ICA-RUS) project funded by the Environment Research and Technology Development Fund (S-10) of the Ministry of the Environment of Japan. Model simulations were performed on the SGI UV20 at the National Institute for Environmental Studies. We gratefully acknowledge the helpful discussions with Kaoru Tachiiri, Tomohiro Hajima, Takashi Arakawa, Junichi Tsutsui, and Michio Kawamiya. The authors are much indebted to Keita Matsumoto, Kuniyasu Hamada, Kenryou Kataumi, Eiichi Hirohashi, Futoshi Takeuchi, Nobuaki Morita, and Kenji Yoshimura at NEC Corporation for their support in model development.

Financial support

This research has been supported by the Ministry of Education, Culture, Sports, Science and Technology of Japan (Integrated Research Program for Advancing Climate Models) (grant no. JPMXD0717935715) and the Ministry of the Environment of Japan (The Environment Research and Technology Development Fund (S-10, grant no. JPMEERF12S11000)).

Review statement

This paper was edited by Olivier Marti and reviewed by two anonymous referees.


Ainsworth, E. A. and Long, P.: What have we learned from 15 years of free-air CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2, New Phytol., 165, 351–372, 2005. 

Alexander, P., Rounsevell, M. D. A., Dislich, C., Dodson, J. R., Engström, K., and Moran, D.: Drivers for global agricultural land use change: The nexus of diet, population, yield and bioenergy, Global Environ. Chang., 35, 138–147, 2015. 

Alexander, P., Prestele, R., Verburg, P. H., Arneth, A., Baranzelli, C., Batista e Silva, F., Brown, C., Butler, A., Calvin, K., Dendoncker, N., Doelman, J. C., Dunford, R., Engström, K., Eitelberg, D., Fujimori, S., Harrison, P. A., Hasegawa, T., Havlik, P., Holzhauer, S., Humpenöder, F., Jacobs-Crisioni, C., Jain, A. K., Krisztin, T., Kyle, P., Lavalle, C., Lenton, T., Liu, J., Meiyappan, P., Popp, A., Powell, T., Sands, R. D., Schaldach, R., Stehfest, E., Steinbuks, J., Tabeau, A., van Meijl, H., Wise, M. A., and Rounsevell, M. D. A.: Assessing uncertainties in land cover projections, Global Change Biol., 23, 767–781, 2017. 

Alexander, P., Rabin, S., Anthoni, P., Henry, R., Pugh, T. A. M., Rounsevell, M. D. A., and Arneth, A.: Adaptation of global land use and management intensity to changes in climate and atmospheric carbon dioxide, Global Change Biol., 24, 2791–2809, 2018. 

Arent, D. J., Döll, P., Strzepek, K. M., Jiménez Cisneros, B. E., Reisinger, A., Tóth, F. L., and Oki, T.: Cross-chapter box on the water–energy–food/feed/fiber nexus as linked to climate change, in: Climate Change 2014: Impacts, Adaptation, and Vulnerability – Part A: Global and Sectoral Aspects, Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel of Climate Change, edited by: Field, C. B., Barros, V. R., Dokken, D. J., Mach, K. J., Mastrandrea, M. D., Bilir, T. E., Chatterjee, M., Ebi, K. L., Estrada, Y. O., Genova, R. C., Girma, B., Kissel, E. S., Levy, A. N., MacCracken, S., Mastrandrea, P. R., and White, L. L., Cambridge University Press, Cambridge, UK, New York, NY, USA, 2014. 

Baldocchi, D.: An analytical solution for coupled leaf photosynthesis and stomatal conductance models, Tree Physiol., 14, 1069–1079, 1994. 

Beer, C., Reichstein, M., Tomelleri, E., Ciais, P., Jung, M., Carvalhais, N., Rödenbeck, C., Arain, M. A., Baldocchi, D., Bonan, G. B., Bondeau, A., Cescatti, A., Lasslop, G., Lindroth, A., Lomas, M., Luyssaert, S., Margolis, H., Oleson, K. W., Roupsard, O., Veenendaal, E., Viovy, N., Williams, C., Woodward, F. I., and Papale, D.: Terrestrial gross carbon dioxide uptake: Global distribution and covariation with climate, Science, 329, 834–838,, 2010. 

Bentsen, M., Bethke, I., Debernard, J. B., Iversen, T., Kirkevåg, A., Seland, Ø., Drange, H., Roelandt, C., Seierstad, I. A., Hoose, C., and Kristjánsson, J. E.: The Norwegian Earth System Model, NorESM1-M – Part 1: Description and basic evaluation of the physical climate, Geosci. Model Dev., 6, 687–720,, 2013. 

Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology / Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant, Hydrol. Sci. Bull., 24, 43–69, 1979. 

BirdLife International: World database of key biodiversity areas. Developed by the KBA partnership, BirdLife International, International Union for the Conservation of Nature, Amphibian Survival Alliance, Conservation International, Critical Ecosystem Partnership Fund, Global Environment Facility, Global Wildlife Conservation, NatureServe, Rainforest Trust, Royal Society for the Protection of Birds, Wildlife Conservation Society and World Wildlife Fund, available at: (last access: 8 January 2018), 2017. 

Bondeau, A., Smith, P. C., Zaehle, S., Schaphoff, S., Lucht, W., Cramer, W., Gerten, D., Lotze-Campen, H., MÜLler, C., Reichstein, M., and Smith, B.: Modelling the role of agriculture for the 20th century global terrestrial carbon balance, Glob. Change Biol., 13, 679–706, 2007. 

Brovkin, V., Boysen, L., Arora, V. K., Boisier, J. P., Cadule, P., Chini, L., Claussen, M., Friedlingstein, P., Gayler, V., van den Hurk, B. J. J. M., Hurtt, G. C., Jones, C. D., Kato, E., de Noblet-Ducoudré, N., Pacifico, F., Pongratz, J., and Weiss, M.: Effect of Anthropogenic Land-Use and Land-Cover Changes on Climate and Land Carbon Storage in CMIP5 Projections for the Twenty-First Century, J. Climate, 26, 6859–6881, 2013. 

Butler, J. H.: Economic Geography: Spatial and Environmental Aspects of Economic Activity, John Wiley, New York, 402 pp., 1980 

Calvin, K. and Bond-Lamberty, B.: Integrated human-earth system modeling–state of the science and future directions, Environ. Res. Lett., 13, 063006,, 2018. 

Chaudhari, S., Pokhrel, Y., Moran, E., and Miguez-Macho, G.: Multi-decadal hydrologic change and variability in the Amazon River basin: understanding terrestrial water storage variations and drought characteristics, Hydrol. Earth Syst. Sci., 23, 2841–2862,, 2019. 

Chen, L. and Dirmeyer, P. A.: Adapting observationally based metrics of biogeophysical feedbacks from land cover/land use change to climate modeling, Environ. Res. Lett., 11, 034002,, 2016. 

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., et al.: Carbon and other biogeochemical cycles, in: Climate Change 2013: The Physical Science Basis Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Stocker, T. F., et al., 465–570, Cambridge, United Kingdom and New York, NY, USA, Cambridge Univ. Press, 2013. 

Collatz, G. J., Ball, J. T., Grivet, C., and Berry, J. A.: Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agr. Forest Meteorol., 54, 107–136, 1991 

Collins, W. D., Craig, A. P., Truesdale, J. E., Di Vittorio, A. V., Jones, A. D., Bond-Lamberty, B., Calvin, K. V., Edmonds, J. A., Kim, S. H., Thomson, A. M., Patel, P., Zhou, Y., Mao, J., Shi, X., Thornton, P. E., Chini, L. P., and Hurtt, G. C.: The integrated Earth system model version 1: formulation and functionality, Geosci. Model Dev., 8, 2203–2219,, 2015. 

Community Earth System Model Project, available at:, last access: 1 July 2019. 

Debele, B., Srinivasan, R., and Parlange, J. Y.: Accuracy evaluation of weather data generation and disaggregation methods at finer timescales, Adv. Water Resour., 30, 1286–1300,, 2007. 

Dietrich, J. P., Bodirsky, B. L., Humpenöder, F., Weindl, I., Stevanović, M., Karstens, K., Kreidenweis, U., Wang, X., Mishra, A., Klein, D., Ambrósio, G., Araujo, E., Yalew, A. W., Baumstark, L., Wirth, S., Giannousakis, A., Beier, F., Chen, D. M.-C., Lotze-Campen, H., and Popp, A.: MAgPIE 4 – a modular open-source framework for modeling global land systems, Geosci. Model Dev., 12, 1299–1317,, 2019. 

Dietrich, J. P., Popp, A., and Lotze-Campen, H.: Reducing the loss of information and gaining accuracy with clustering methods in a global land-use model, Ecol. Model., 263, 233–243, 2013. 

Dufresne, J.-L., Foujols, M.-A., Denvil, S., Caubel, A., Marti, O., Aumont, O., Balkanski, Y., Bekki, S., Bellenger, H., Benshila, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., de Noblet, N., Duvel, J.-P., Ethé, C., Fairhead, L., Fichefet, T., Flavoni, S., Friedlingstein, P., Grandpeix, J.-Y., Guez, L., Guilyardi, E., Hauglustaine, D., Hourdin, F., Idelkadi, A., Ghattas, J., Joussaume, S., Kageyama, M., Krinner, G., Labetoulle, S., Lahellec, A., Lefebvre, M.-P., Lefevre, F., Levy, C., Li, Z. X., Lloyd, J., Lott, F., Madec, G., Mancip, M., Marchand, M., Masson, S., Meurdesoif, Y., Mignot, J., Musat, I., Parouty, S., Polcher, J., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Talandier, C., Terray, P., Viovy, N., and Vuichard, N.: Climate change projections using the IPSL-CM5 Earth System Model: from CMIP3 to CMIP5, Clim. Dynam., 40, 2123–2165, 2013. 

Dunne, J. P., John, J. G., Adcroft, A. J., Griffies, S. M., Hallberg, R. W., Shevliakova, E., Stouffer, R. J., Cooke, W., Dunne, K. A., Harrison, M. J., Krasting, J. P., Malyshev, S. L., Milly, P. C. D., Phillipps, P. J., Sentman, L. T., Samuels, B. L., Spelman, M. J., Winton, M., Wittenberg, A. T., and Zadeh, N.: GFDL's ESM2 Global Coupled Climate–Carbon Earth System Models. Part I: Physical Formulation and Baseline Simulation Characteristics, J. Climate, 25, 6646–6665, 2012. 

Ejiri, Y.: A consideration of the Comparative Cost Model Using Three-Dimensional Diagrams, Forest Resource Management and Mathematical Modeling, FORMATH, 7, 135–159, 2008. 

Engström, K., Olin, S., Rounsevell, M. D. A., Brogaard, S., van Vuuren, D. P., Alexander, P., Murray-Rust, D., and Arneth, A.: Assessing uncertainties in global cropland futures using a conditional probabilistic modelling framework, Earth Syst. Dynam., 7, 893–915,, 2016. 

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. 

Famien, A. M., Janicot, S., Ochou, A. D., Vrac, M., Defrance, D., Sultan, B., and Noël, T.: A bias-corrected CMIP5 dataset for Africa using the CDF-t method – a contribution to agricultural impact studies, Earth Syst. Dynam., 9, 313–338,, 2018. 

Farquhar, G. D., von Caemmerer, S., and Berry, J. A.: A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species, Planta, 149, 78–90, 1980. 

Feddema, J. J., Oleson, K. W., Bonan, G. B., Mearns, L. O., Buja, L. E., Meehl, G. A., and Washington, W. M.: The Importance of Land-Cover Change in Simulating Future Climates, Science, 310, 1674–1678, 2005. 

Felfelani, F., Wada, Y., Longuevergne, L., and Pokhrel, Y.: Natural and human-induced terrestrial water storage change: A global analysis using hydrological models and GRACE, J. Hydrol., 553, 105–118, 2017. 

Findell, K. L., Berg, A., Gentine, P., Krasting, J. P., Lintner, B. R., Malyshev, S., Santanello, J. A., and Shevliakova, E.: The impact of anthropogenic land use and land cover change on regional climate extremes, Nat. Commun., 8, 989,, 2017. 

Foley, J. A., Ramankutty, N., Brauman, K. A., Cassidy, E. S., Gerber, J. S., Johnston, M., Mueller, N. D., O'Connell, C., Ray, D. K., West, P. C., Balzer, C., Bennett, E. M., Carpenter, S. R., Hill, J., Monfreda, C., Polasky, S., Rockström, J., Sheehan, J., Siebert, S., Tilman, D., and Zaks, D. P. M.: Solutions for a cultivated planet, Nature, 478, 337–342, 2011. 

Food and Agriculture Organization of the United Nations (FAO): FAOSTAT, Rome, Italy: The Stat. Div. of FAO, available at: (last access: 20 July 2020), 2019. 

Friedl, M. A., Sulla-Menashe, D., Tan, B., Schneider, A., Ramankutty, N., Sibley, A., and Huang, X.: MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets, Remote Sens. Environ., 114, 168–182, 2010. 

Fujimori, S., Masui, T., Matsuoka, Y.: AIM/CGE [basic] manual, Discussion Paper Series, 1–74, 2012. 

Fujimori, S., Abe, M., Kinoshita, T., Hasegawa, T., Kawase, H., Kushida, K., Masui, T., Oka, K., Shiogama, H., Takahashi, K., Tatebe, H., and Yoshikawa, M.: Downscaling Global Emissions and Its Implications Derived from Climate Model Experiments, PLOS ONE 12, e0169733,, 2017a. 

Fujimori, S., Hasegawa, T., Masui, T., Takahashi, K., Herran, D. S., Dai, H., Hijioka, Y., and Kainuma, M.: SSP3: AIM implementation of Shared Socioeconomic Pathways, Global Environ. Chang., 42, 268–283,, 2017b. 

Goudriaan, J. and van Laar, H. H.: Modelling potential crop growth processes, Kluwer Academic Publishers, 1994. 

Hajima, T., Watanabe, M., Yamamoto, A., Tatebe, H., Noguchi, M. A., Abe, M., Ohgaito, R., Ito, A., Yamazaki, D., Okajima, H., Ito, A., Takata, K., Ogochi, K., Watanabe, S., and Kawamiya, M.: Development of the MIROC-ES2L Earth system model and the evaluation of biogeochemical processes and feedbacks, Geosci. Model Dev., 13, 2197–2244,, 2020. 

Hanasaki, N., Kanae, S., and Oki, T.: A reservoir operation scheme for global river routing models, J. Hydrol., 327, 22–41, 2006. 

Hanasaki, N., Kanae, S., Oki, T., Masuda, K., Motoya, K., Shirakawa, N., Shen, Y., and Tanaka, K.: An integrated model for the assessment of global water resources – Part 1: Model description and input meteorological forcing, Hydrol. Earth Syst. Sci., 12, 1007–1025,, 2008a. 

Hanasaki, N., Kanae, S., Oki, T., Masuda, K., Motoya, K., Shirakawa, N., Shen, Y., and Tanaka, K.: An integrated model for the assessment of global water resources – Part 2: Applications and assessments, Hydrol. Earth Syst. Sci., 12, 1027–1037,, 2008b. 

Hanasaki, N., Inuzuka, T., Kanae, S., and Oki, T.: An estimation of global virtual water flow and sources of water withdrawal for major crops and livestock products using a global hydrological model, J. Hydrol., 384, 232–244, 2010. 

Hasegawa, T., Fujimori, S., Ito, A., Takahashi, K., and Masui, T.: Global land-use allocation model linked to an integrated assessment model, Sci. Total Environ., 580, 787–796, 2017. 

Havlík, P., Schneider, U. A., Schmid, E., Böttcher, H., Fritz, S., Skalský, R., Aoki, K., Cara, S. D., Kindermann, G., Kraxner, F., Leduc, S., McCallum, I., Mosnier, A., Sauer, T., and Obersteiner, M.: Global land-use implications of first and second generation biofuel targets, Energy Policy, 39, 5690–5702, 2011. 

Hempel, S., Frieler, K., Warszawski, L., Schewe, J., and Piontek, F.: A trend-preserving bias correction – the ISI-MIP approach, Earth Syst. Dynam., 4, 219–236,, 2013. 

Hibbard, K., Janetos, A., van Vuuren, D. P., Pongratz, J., Rose, S. K., Betts, R., Herold, M., and Feddema, J. J.: Research priorities in land use and land-cover change for the Earth system and integrated assessment modelling, Int. J. Climatol., 30, 2118–2128, 2010. 

Hirsch, A. L., Guillod, B. P., Seneviratne, S. I., Beyerle, U., Boysen, L. R., Brovkin, V., Davin, E. L., Doelman, J. C., Kim, H., Mitchell, D. M., Nitta, T., Shiogama, H., Sparrow, S., Stehfest, E., Vuuren, D. P., and Wilson, S.: Biogeophysical Impacts of Land-Use Change on Climate Extremes in Low-Emission Scenarios: Results From HAPPI-Land, Earth's Future, 6, 396–409, 2018. 

Howden, S. M., Soussana, J.-F., Tubiello, F. N., Chhetri, N., Dunlop, M., and Meinke, H.: Adapting agriculture to climate change, P. Natl. Acad. Sci. USA, 104, 19691–19696, 2007. 

Humpenöder, F., Popp, A., Stevanovic, M., Müller, C., Bodirsky, B. L., Bonsch, M., Dietrich, J. P., Lotze-Campen, H., Weindl, I., Biewald, A., and Rolinski, S.: Land-Use and Carbon Cycle Responses to Moderate Climate Change: Implications for Land-Based Mitigation?, Environ. Sci. Tech., 49, 6731–6739, 2015. 

Hurtt, G. C., Frolking, S., Fearon, M. G., Moore, B., Shevliakova, E., Malyshev, S., Pacala, S. W., and Houghton, R. A.: The underpinnings of land-use history: three centuries of global gridded land-use transitions, wood-harvest activity, and resulting secondary lands, Glob. Change Biol., 12, 1208–1229, 2006. 

Iizumi, T., Yokozawa, M., Sakurai, G., Travasso, M.I., Romanernkov, V., Oettli, P., Newby, T., Ishigooka, Y., and Furuya, J.: Historical changes in global yields: Major cereal and legume crops from 1982 to 2006, Global Ecol. Biogeogr., 23, 346–357, 2013. 

Iizumi, T., Okada, M., and Yokozawza, M.: A meteorological forcing data set for global crop modeling: Development, evaluation, and intercomparison, J. Geophys. Res.-Atmos., 119, 363–384,, 2014. 

Iizumi, T.: GDHY, Data set, Data Integration and Analysis System (DIAS),, 2017. 

Ito, A.: Changing ecophysiological processes and carbon budget in East Asian ecosystems under near-future changes in climate: Implications for long-term monitoring from a process-based model, J. Plant Res., 123, 577–588, 2010. 

Ito, A. and Inatomi, M.: Use of a process-based model for assessing the methane budgets of global terrestrial ecosystems and evaluation of uncertainty, Biogeosciences, 9, 759–773,, 2012. 

Ito, A., Nishina, K., Ishijima, K., Hashimoto, S., and Inatomi, M.: Emissions of nitrous oxide (N2O) from soil surfaces and their historical changes in East Asia: a model-based assessment, Progr. Earth Planet. Sci., 5, 55,, 2018. 

Ito, A., Nishina, K., Reyer, C. P. O., François, L., Henrot, A.-J., Munhoven, G., Jacquemin, I., Tian, H., Yang, J., Pan, S., Morfopoulos, C., Betts, R., Hickler, T., Steinkamp, J., Ostberg, S., Schaphoff, S., Ciais, P., Chang, J., Rafique, R., Zeng, F., and Zhao, F.: Photosynthetic productivity and its efficiencies in ISIMIP2a biome models: benchmarking for impact assessment studies, Environ. Res. Lett., 12, 085001,, 2017. 

IUCN and UNEP‐WCMC: The World Database on Protected Areas (WDPA), UNEP‐WCMC, Cambridge, UK, available at: (last access: 2 October 2020), 2018. 

Jones, C. D., Hughes, J. K., Bellouin, N., Hardiman, S. C., Jones, G. S., Knight, J., Liddicoat, S., O'Connor, F. M., Andres, R. J., Bell, C., Boo, K.-O., Bozzo, A., Butchart, N., Cadule, P., Corbin, K. D., Doutriaux-Boucher, M., Friedlingstein, P., Gornall, J., Gray, L., Halloran, P. R., Hurtt, G., Ingram, W. J., Lamarque, J.-F., Law, R. M., Meinshausen, M., Osprey, S., Palin, E. J., Parsons Chini, L., Raddatz, T., Sanderson, M. G., Sellar, A. A., Schurer, A., Valdes, P., Wood, N., Woodward, S., Yoshioka, M., and Zerroukat, M.: The HadGEM2-ES implementation of CMIP5 centennial simulations, Geosci. Model Dev., 4, 543–570,, 2011. 

K-1 model developers: K-1 coupled GCM (MIROC) description, K-1 Tech. Rep., 1, edited by: Hasumi, H. and Emori, S., Center for Climate System Research, the University of Tokyo, Tokyo, 34 pp., 2004. 

Kaschuk, G., Yin, X., Hungria, M., Leffelaar, P. A., Giller, K. E., and Kuyper, T. W.: Photosynthetic adaptation of soybean due to varying effectiveness of N2 fixation by two distinct Bradyrhizobium japonicum strains, Environ. Exp. Bot., 76, 1–6, 2012. 

Kato, E. and Yamagata, Y.: BECCS capability of dedicated bioenergy crops under a future land-use scenario targeting net negative carbon emissions, Earth's Future, 2, 421–439, 2014. 

Koirala, S., Yeh, P. J. F., Hirabayashi, Y., Kanae, S., and Oki, T.: Global-scale land surface hydrologic modeling with the representation of water table dynamics, J. Geophys. Res.-Atmos., 119, 75–89, 2014. 

Kriegler, E. and Lucht, W.: Overview of the PIK REMINDMAgPIE-LPJml integrated assessment framework, available at: (last access: 1 July 2019), 2015. 

Krysanova, V., Müller-Wohlfeil, D.-I., and Becker, A.: Development and test of a spatially distributed hydrological/water quality model for mesoscale watersheds, Ecol. Model., 106, 261–289, 1998. 

Lawrence, D. M., Hurtt, G. C., Arneth, A., Brovkin, V., Calvin, K. V., Jones, A. D., Jones, C. D., Lawrence, P. J., de Noblet-Ducoudré, N., Pongratz, J., Seneviratne, S. I., and Shevliakova, E.: The Land Use Model Intercomparison Project (LUMIP) contribution to CMIP6: rationale and experimental design, Geosci. Model Dev., 9, 2973–2998,, 2016. 

Le Quéré, C., Andrew, R. M., Friedlingstein, P., Sitch, S., Hauck, J., Pongratz, J., Pickers, P. A., Korsbakken, J. I., Peters, G. P., Canadell, J. G., Arneth, A., Arora, V. K., Barbero, L., Bastos, A., Bopp, L., Chevallier, F., Chini, L. P., Ciais, P., Doney, S. C., Gkritzalis, T., Goll, D. S., Harris, I., Haverd, V., Hoffman, F. M., Hoppema, M., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Johannessen, T., Jones, C. D., Kato, E., Keeling, R. F., Goldewijk, K. K., Landschützer, P., Lefèvre, N., Lienert, S., Liu, Z., Lombardozzi, D., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S., Neill, C., Olsen, A., Ono, T., Patra, P., Peregon, A., Peters, W., Peylin, P., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rocher, M., Rödenbeck, C., Schuster, U., Schwinger, J., Séférian, R., Skjelvan, I., Steinhoff, T., Sutton, A., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Laan-Luijkx, I. T., van der Werf, G. R., Viovy, N., Walker, A. P., Wiltshire, A. J., Wright, R., Zaehle, S., and Zheng, B.: Global Carbon Budget 2018, Earth Syst. Sci. Data, 10, 2141–2194,, 2018. 

Lehner, B., Liermann, C. R., Revenga, C., Vörösmarty, C., Fekete, B., Crouzet, P., Döll, P., Endejan, M., Frenken, K., Magome, J., Nilsson, C., Robertson, J. C., Rödel, R., Sindorf, N., and Wisser, D.: High-resolution mapping of the world's reservoirs and dams for sustainable river-flow management, Front. Ecol. Environ., 9, 494–502, 2011. 

Letourneau, A., Verburg, P. H., and Stehfest, E.: A land-use systems approach to represent land-use dynamics at continental and global scales, Environ. Model. Softw., 33, 61–79, 2012. 

Liu, B., Asseng, S., Müller, C., Ewert, F., Elliott, J., Lobell, David B., Martre, P., Ruane, Alex C., Wallach, D., Jones, James W., Rosenzweig, C., Aggarwal, Pramod K., Alderman, Phillip D., Anothai, J., Basso, B., Biernath, C., Cammarano, D., Challinor, A., Deryng, D., Sanctis, Giacomo D., Doltra, J., Fereres, E., Folberth, C., Garcia-Vila, M., Gayler, S., Hoogenboom, G., Hunt, Leslie A., Izaurralde, Roberto C., Jabloun, M., Jones, Curtis D., Kersebaum, Kurt C., Kimball, Bruce A., Koehler, A.-K., Kumar, Soora N., Nendel, C., O'Leary, Garry J., Olesen, Jørgen E., Ottman, Michael J., Palosuo, T., Prasad, P. V. V., Priesack, E., Pugh, Thomas A. M., Reynolds, M., Rezaei, Ehsan E., Rötter, Reimund P., Schmid, E., Semenov, Mikhail A., Shcherbak, I., Stehfest, E., Stöckle, Claudio O., Stratonovitch, P., Streck, T., Supit, I., Tao, F., Thorburn, P., Waha, K., Wall, Gerard W., Wang, E., White, Jeffrey W., Wolf, J., Zhao, Z., and Zhu, Y.: Similar estimates of temperature impacts on global wheat yield by three independent methods, Nat. Clim. Change, 6, 1130–1136, 2016. 

Lotze-Campen, H., Müller, C., Bondeau, A., Rost, S., Popp, A., and Lucht, W.: Global food demand, productivity growth, and the scarcity of land and water resources: a spatially explicit mathematical programming approach, Agr. Econom., 39, 325–338,, 2008. 

Mahmood, R., Pielke Sr., R. A., Hubbard, K. G., Niyogi, D., Dirmeyer, P. A., McAlpine, C., Carleton, A. M., Hale, R., Gameda, S., Beltrán-Przekurat, A., Baker, B., McNider, R., Legates, D. R., Shepherd, M., Du, J., Blanken, P. D., Frauenfeld, O. W., Nair, U. S., and Fall, S.: Land cover changes and their biogeophysical effects on climate, Int. J. Climatol., 34, 929–953, 2014. 

McGuire, A., Sitch, D. S., Clein, J. S., Dargaville, R., Esser, G., Foley, J., Heimann, M., Joos, F., Kaplan, J., Kicklighter, D. W., Meier, R. A., Melillo, J. M., Moore III, B., Prentice, I. C., Ramankutty, N., Reichenau, T., Schloss, A., Tian, H., Williams, L. J., and Wittenberg, U.: Carbon balance of the terrestrial biosphere in the twentieth century: analysis of CO2, climate and land use effects with four process-based ecosystem models, Glob. Biogeochem. Cy., 15, 183–206, 2001. 

Medlyn, B. E., Dreyer, E., Ellsworth, D., Forstreuter, M., Harley, P. C., Kirschbaum, M. U. F., Le Roux, X., Montpied, P., Strassemeyer, J., Walcroft, A., Wang, K., and Loustau, D.: Temperature response of parameters of a biochemically based model of photosynthesis. II. A review of experimental data, Plant Cell Environ., 25, 1167–1179, 2002. 

Meiyappan, P., Dalton, M., O'Neill, B. C., and Jain, A. K.: Spatial modeling of agricultural land use change at global scale, Ecol. Model., 291, 152–174, 2014. 

Monfreda, C., Ramankutty, N., and Foley, J.: Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000, Global Biogeochem. Cy., 22, GB1022,, 2008. 

Monsi, M. and Saeki, T.: Über den Lichtfaktor in den Pflanzengesellschaften und seine Bedeutung für die Stoffproduktion, Jpn. J. Bot., 14, 22–52, 1953. 

Mori, S., Washida, T., Kurosawa, A., and Masui, T.: Assessment of mitigation strategies as tools for risk management under future uncertainties: a multi-model approach, Sustain. Sci., 13, 329–349,, 2018. 

Moss, R. H., Edmonds, J. A., Hibbard, K. A., Manning, M. R., Rose, S. K., van Vuuren, D. P., Carter, T. R., Emori, S., Kainuma, M., Kram, T., Meehl, G. A., Mitchell, J. F. B., Nakicenovic, N., Riahi, K., Smith, S. J., Stouffer, R. J., Thomson, A. M., Weyant, J. P., and Wilbanks, T. J.: The next generation of scenarios for climate change research and assessment, Nature, 463, 747–756, 2010. 

Müller, C., Elliott, J., Chryssanthacopoulos, J., Arneth, A., Balkovic, J., Ciais, P., Deryng, D., Folberth, C., Glotter, M., Hoek, S., Iizumi, T., Izaurralde, R. C., Jones, C., Khabarov, N., Lawrence, P., Liu, W., Olin, S., Pugh, T. A. M., Ray, D. K., Reddy, A., Rosenzweig, C., Ruane, A. C., Sakurai, G., Schmid, E., Skalsky, R., Song, C. X., Wang, X., de Wit, A., and Yang, H.: Global gridded crop model evaluation: benchmarking, skills, deficiencies and implications, Geosci. Model Dev., 10, 1403–1422,, 2017. 

Müller-Hansen, F., Schlüter, M., Mäs, M., Donges, J. F., Kolb, J. J., Thonicke, K., and Heitzig, J.: Towards representing human behavior and decision making in Earth system models – an overview of techniques and approaches, Earth Syst. Dynam., 8, 977–1007,, 2017. 

Murakami, D. and Yamagata, Y.: Estimation of Gridded Population and GDP Scenarios with Spatially Explicit Statistical Downscaling, Sustainability, 11, 2106,, 2019. 

Neitsch, S. L., Arnold, J. G., Kiniry, J. R., Williams, J. R., and King, K. W.: Soil and water assessment tool theoretical documentation (version 2005), United States Department of Agriculture, 2005. 

Nitta, T., Yoshimura, K., Takata, K., O'ishi, R., Sueyoshi, T., Kanae, S., Oki, T., Abe-Ouchi, A., and Liston, G. E.: Representing Variability in Subgrid Snow Cover and Snow Depth in a Global Land Model: Offline Validation, J. Climate, 27, 3318–3330, 2014. 

O'Neill, B. C., Kriegler, E., Ebi, K. L., Kemp-Benedict, E., Riahi, K., Rothman, D. S., van Ruijven, B. J., van Vuuren, D. P., Birkmann, J., Kok, K., Levy, M., and Solecki, W.: The roads ahead: Narratives for shared socioeconomic pathways describing world futures in the 21st century, Global Environ. Chang., 42, 169–180, 2017. 

Okada, M., Iizumi, T., Sakurai, G., Hanasaki, N., Sakai, T., Okamoto, K., and Yokozawa, M.: Modeling irrigation-based climate change adaptation in agriculture: Model development and evaluation in Northeast China, J. Adv. Model. Earth Syst., 7, 1409–1424, 2015 

Oki, T. and Sud, Y. C.: Design of Total Runoff Integrating Pathways (TRIP) – A Global River Channel Network, Earth Interact., 2, 1–37, 1998. 

Olin, S., Schurgers, G., Lindeskog, M., Wårlind, D., Smith, B., Bodin, P., Holmér, J., and Arneth, A.: Modelling the response of yields and tissue C : N to changes in atmospheric CO2 and N management in the main wheat regions of western Europe, Biogeosciences, 12, 2489–2515,, 2015. 

Parry, M. L., Rosenzweig, C., Iglesias, A., Livermore, M., and Fischer, G.: Effects of climate change on global food production under SRES emissions and socio-economic scenarios, Global Environ. Chang., 14, 53–67, 2004. 

Pokhrel, Y., Hanasaki, N., Wada, Y., and Kim, H.: Recent progresses in incorporating human land–water management into global land surface models toward their integration into Earth system models, WIREs Water, 3, 548–574, 2016. 

Pokhrel, Y., Felfelani, F., Shin, S., Yamada, T. J., and Satoh, Y.: Modeling large-scale human alteration of land surface hydrology and climate, Geosci. Lett., 4, 10,, 2017. 

Pokhrel, Y., Hanasaki, N., Yeh, P. J. F., Yamada, T. J., Kanae, S., and Oki, T.: Model estimates of sea-level change due to anthropogenic impacts on terrestrial water storage, Nat. Geosci., 5, 389–392, 2012a. 

Pokhrel, Y., Koirala, S., Yeh, P. J. F., Hanasaki, N., Longuevergne, L., Kanae, S., and Oki, T.: Incorporation of groundwater pumping in a global Land Surface Model with the representation of human impacts, Water Resour. Res., 51, 78–96, 2015. 

Pokhrel, Y., Hanasaki, N., Koirala, S., Cho, J., Yeh, P. J. F., Kim, H., Kanae, S., and Oki, T.: Incorporating Anthropogenic Water Regulation Modules into a Land Surface Model, J. Hydrometeorol., 13, 255–269, 2012b. 

Popp, A., Dietrich, J. P., Lotze-Campen, H., Klein, D., Bauer, N., Krause, M., Beringer, T., Gerten, D., and Edenhofer, O.: The economic potential of bioenergy for climate change mitigation with special attention given to implications for the land system, Environ. Res. Lett., 6, 034017,, 2011. 

Popp, A., Calvin, K., Fujimori, S., Havlik, P., Humpenöder, F., Stehfest, E., Bodirsky, B. L., Dietrich, J. P., Doelmann, J. C., Gusti, M., Hasegawa, T., Kyle, P., Obersteiner, M., Tabeau, A., Takahashi, K., Valin, H., Waldhoff, S., Weindl, I., Wise, M., Kriegler, E., Lotze-Campen, H., Fricko, O., Riahi, K., and Vuuren, D. P. v.: Land-use futures in the shared socio-economic pathways, Global Environ. Chang., 42, 331–345, 2017. 

Porter, J. R., Xie, L., Challinor, A. J., Cochrane, K., Howden, S. M., Iqbal, M. M., Lobell, D. B., and Travasso, M. I.: Food security and food production systems, in: Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects, Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel of Climate Change, edited by: Field, C. B., Barros, V. R., Dokken, D. J., Mach, K. J., Mastrandrea, M. D., Bilir, T. E., Chatterjee, M., Ebi, K. L., Estrada, Y. O., Genova, R. C., Girma, B., Kissel, E. S., Levy, A. N., MacCracken, S., Mastrandrea, P. R., and White, L. L., Cambridge University Press, Cambridge, UK, New York, NY, USA, 2014. 

Pugh, T. A. M., Müller, C., Elliott, J., Deryng, D., Folberth, C., Olin, S., Schmid, E., and Arneth, A.: Climate analogues suggest limited potential for intensification of production on current croplands under climate change, Nat. Commun., 7, 12608,, 2016. 

Richards, L. A.: Capillary Conduction Of Liquids Through Porous Mediums, J. Appl. Phys., 1, 318–333, 1931. 

Robinson, D. T., Di Vittorio, A., Alexander, P., Arneth, A., Barton, C. M., Brown, D. G., Kettner, A., Lemmen, C., O'Neill, B. C., Janssen, M., Pugh, T. A. M., Rabin, S. S., Rounsevell, M., Syvitski, J. P., Ullah, I., and Verburg, P. H.: Modelling feedbacks between human and natural processes in the land system, Earth Syst. Dynam., 9, 895–914,, 2018. 

Rolinski, S., Müller, C., Heinke, J., Weindl, I., Biewald, A., Bodirsky, B. L., Bondeau, A., Boons-Prins, E. R., Bouwman, A. F., Leffelaar, P. A., te Roller, J. A., Schaphoff, S., and Thonicke, K.: Modeling vegetation and carbon dynamics of managed grasslands at the global scale with LPJmL 3.6, Geosci. Model Dev., 11, 429–451,, 2018. 

Romero-Lankao, P., Smith, J. B., Davidson, D. J., Diffenbaugh, N. S., Kinney, P. L., Kirshen, P., Kovacs, P., and Villers-Ruiz, L.: North America, in: Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part B: Regional Aspects, Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel of Climate Change, edited by: Barros, V. R., Field, C. B., Dokken, D. J., Mastrandrea, M. D., Mach, K. J., Bilir, T. E., Chatterjee, M., Ebi, K. L., Estrada, Y. O., Genova, R. C., Girma, B., Kissel, E. S., Levy, A. N., MacCracken, S., Mastrandrea, P. R., and White, L. L., Cambridge University Press, Cambridge, UK, New York, NY, USA, 2014. 

Rosenzweig, C., Elliott, J., Deryng, D., Ruane, A. C., Müller, C., Arneth, A., Boote, K. J., Folberth, C., Glotter, M., Khabarov, N., Neumann, K., Piontek, F., Pugh, T. A. M., Schmid, E., Stehfest, E., Yang, H., and Jones, J. W.: Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison, P. Natl. Acad. Sci. USA, 111, 3268–3273, 2014. 

Rounsevell, M. D. A., Arneth, A., Alexander, P., Brown, D. G., de Noblet-Ducoudré, N., Ellis, E., Finnigan, J., Galvin, K., Grigg, N., Harman, I., Lennox, J., Magliocca, N., Parker, D., O'Neill, B. C., Verburg, P. H., and Young, O.: Towards decision-based global land use models for improved understanding of the Earth system, Earth Syst. Dynam., 5, 117–137,, 2014. 

Sacks, W. J., Deryng, D., Foley, J. A., and Ramankutty, N.: Crop planting dates: an analysis of global patterns, Global Ecol. Biogeogr., 19, 607–620, 2010. 

Sakurai, G., Iizumi, T., Nishimori, M., and Yokozawa, M.: How much has the increase in atmospheric CO2 directly affected past soybean production?, Sci. Rep., 4, 4978,, 2014. 

Save, H., Bettadpur, S., and Tapley, B. D.: High-resolution CSR GRACE RL05 mascons, J. Geophys. Res.-Sol. Ea., 121, 7547–7569, 2016. 

Scanlon, B. R., Zhang, Z., Save, H., Sun, A. Y., Schmied, H. M., van Beek, L. P., Wiese, D. N., Wada, Y., Long, D., and Reedy, R. C.: Global models underestimate large decadal declining and rising water storage trends relative to GRACE satellite data, P. Natl. Acad. Sci. USA, 115, E1080–E1089,, 2018. 

Scholes, R. J. and Brown de Colstoun, E.: ISLSCP II Global Gridded Soil Characteristics, in: Hall, F. G., Collatz, G., Meeson, B., Los, S., Brown de Colstoun, E., and Landis, D., ISLSCP Initiative II Collection, Oak Ridge National Laboratory Distributed Active Archive Center, Oak Ridge, Tennessee, U.S.A., available at: (last access: 1 October 2020),, 2011. 

Sellers, P. J., Randall, D. A., Collatz, G. J., Berry, J. A., Field, C. B., Dazlich, D. A., Zhang, C., Collelo, G. D., and Bounoua, L.: A revised land surface parameterization (SiB2) for atmospheric GCMs, Part I: Model Formulation, J. Climate, 9, 676–705, 1996a. 

Sellers, P. J., Tucker, C. J., Collatz, G. J., Los, S. O., Justice, C. O., Dazlich, D. A., and Randall, D. A.: A revised land surface parameterization (SiB2) for atmospheric GCMs, Part II: The generation of global fields of terrestrial biophysical parameters from satellite data, J. Climate, 9, 706–737, 1996b. 

Smith, P., Gregory, P. J., van Vuuren, D., Obersteiner, M., Havlík, P., Rounsevell, M., Woods, J., Stehfest, E., and Bellarby, J.: Competition for land, Philos. T. R. Soc. B, 365, 2941–2957, 2010. 

Smith, P., Haberl, H., Popp, A., Erb, K.-h., Lauk, C., Harper, R., Tubiello, F. N., de Siqueira Pinto, A., Jafari, M., Sohi, S., Masera, O., Böttcher, H., Berndes, G., Bustamante, M., Ahammad, H., Clark, H., Dong, H., Elsiddig, E. A., Mbow, C., Ravindranath, N. H., Rice, C. W., Robledo Abad, C., Romanovskaya, A., Sperling, F., Herrero, M., House, J. I., and Rose, S.: How much land-based greenhouse gas mitigation can be achieved without compromising food security and environmental goals?, Glob. Change Biol., 19, 2285–2302, 2013. 

Stieglitz, M., Rind, D., Famiglietti, J., and Rosenzweig, C.: An Efficient Approach to Modeling the Topographic Control of Surface Hydrology for Regional and Global Climate Modeling, J. Climate, 10, 118–137, 1997. 

Takata, K., Emori, S., and Watanabe, T.: Development of the minimal advanced treatments of surface interaction and runoff, Global Planet. Change, 38, 209–222, 2003 

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498, 2012. 

Thornton, P. E., Calvin, K., Jones, A. D., Di Vittorio, A. V., Bond-Lamberty, B., Chini, L., Shi, X., Mao, J., Collins, W. D., Edmonds, J., Thomson, A., Truesdale, J., Craig, A., Branstetter, M. L., and Hurtt, G.: Biospheric feedback effects in a synchronously coupled model of human and Earth systems, Nat. Clim. Change, 7, 496–500, 2017. 

Tilman, D., Balzer, C., Hill, J., and Befort, B. L.: Global food demand and the sustainable intensification of agriculture, P. Natl. Acad. Sci. USA, 108, 20260–20264,, 2011 

United States Department of Agriculture (USDA): Major world crop areas and climatic profiles, Washington, DC, USDA, p. 279, 1994. 

van Vuuren, D. P., Batlle Bayer, L., Chuwah, C., Ganzeveld, L., Hazeleger, W., van den Hurk, B., van Noije, T., O'Neill, B., and Strengers, B. J.: A comprehensive view on climate change: coupling of earth system and integrated assessment models, Environ. Res. Lett., 7, 024012,, 2012. 

van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.-F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J., and Rose, S. K.: The representative concentration pathways: an overview, Clim. Change, 109, 5–31, 2011. 

von Bloh, W., Schaphoff, S., Müller, C., Rolinski, S., Waha, K., and Zaehle, S.: Implementing the nitrogen cycle into the dynamic global vegetation, hydrology, and crop growth model LPJmL (version 5.0), Geosci. Model Dev., 11, 2789–2812,, 2018. 

Verdin, K. L. and Greenlee, S. K.: Development of continental scale digital elevation models and extraction of hydrographic features, in: Proceedings, Third International Conference/Workshop on Integrating GIS and Environmental Modeling, Santa Fe, New Mexico, 21–26 January, 1996, National Center for Geographic Information and Analysis, Santa Barbara, California, 1996. 

Veldkamp, T. I. E., Zhao, F., Ward, P. J., de Moel, H., Aerts, J. C., Schmied, H. M., Portmann, F. T., Masaki, Y., Pokhrel, Y., and Liu, X.: Human impact parameterizations in global hydrological models improve estimates of monthly discharges and hydrological extremes: a multi-model validation study, Environ. Res. Lett., 13, 055008,, 2018. 

Vrugt, J. A., ter Braak, C. J. F., Diks, C. G. H., Robinson, B. A., Hyman, J. M., and Higdon, D.: Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling, Int. J. Nonlinear Sci., 10, 271–288, 2009. 

Watanabe, M., Suzuki, T., O'ishi, R., Komuro, Y., Watanabe, S., Emori, S., Takemura, T., Chikira, M., Ogura, T., Sekiguchi, M., Takata, K., Yamazaki, D., Yokohata, T., Nozawa, T., Hasumi, H., Tatebe, H., and Kimoto, M.: Improved Climate Simulation by MIROC5: Mean States, Variability, and Climate Sensitivity, J. Climate, 23, 6312–6335, 2010. 

Watanabe, S., Hajima, T., Sudo, K., Nagashima, T., Takemura, T., Okajima, H., Nozawa, T., Kawase, H., Abe, M., Yokohata, T., Ise, T., Sato, H., Kato, E., Takata, K., Emori, S., and Kawamiya, M.: MIROC-ESM 2010: model description and basic results of CMIP5-20c3m experiments, Geosci. Model Dev., 4, 845–872,, 2011. 

Watanabe, T.: Bulk parameterization for a vegetated surface and its application to a simulation of nocturnal drainage flow, Bound.-Lay. Meteorol., 70, 13–35, 1994. 

Watkins, M. M., Wiese, D. N., Yuan, D.-N., Boening, C., and Landerer, F. W.: Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons, J. Geophys. Res.-Sol. Ea., 120, 2648–2671, 2015.  

Weinzettel, J., Hertwich, E. G., Peters, G. P., Steen-Olsen, K., and Galli, A.: Affluence drives the global displacement of land use, Global Environ. Chang., 23, 433–438, 2013. 

Wiese, D. N., Yuan, D.-N., Boening, C., Landerer, F. W., and Watkins, M. M.: JPL GRACE Mascon Ocean, Ice, and Hydrology Equivalent Water Height RL05M.1 CRI Filtered Version 2, PO.DAAC, CA, USA, Dataset,, 2016. 

Willett, K. M., Gillett, N. P., Jones, P. D., and Thorne, P. W.: Attribution of observed surface humidity change to human influence, Nature, 449, 710–712, 2007. 

Wiscombe, W. J. and Warren, S. G.: A model for the spectral albedo of snow. I: Pure snow, J. Atmos. Sci., 37, 2712–2733, 1980. 

Wise, M. and Calvin, K.: GCAM3.0 Agriculture and Land Use: Technical Description of Modeling Approach, PNNL-20971, Pacific Northwest National Laboratory, 2011. 

Wise, M., Calvin, K., Kyle, G., Luckow, P., and Edmonds, J.: Economic and physical modeling of land use in GCAM 3.0 and an application to agricultural productivity, land, and terrestrial carbon, Clim. Change Econom., 5, 1450003,, 2014. 

Wu, W., Hasegawa, T., Ohashi, H., Hanasaki, N., Liu, J., Matsui, T., Fujimori, S., Masui, T., and Takahashi, K.: Global advanced bioenergy potential under environmental protection policies and societal transformation measures, GCB Bioenergy, 11, 1041–1055, 2019. 

Yokohata, T., Tanaka, K., Nishina, K., Takahashi, K., Emori, S., Kiguchi, M., Iseri, Y., Honda, Y., Okada, M., Masaki, Y., Yamamoto, A., Shigemitsu, M., Yoshimori, M., Sueyoshi, T., Iwase, K., Hanasaki, N., Ito, A., Sakurai, G., Iizumi, T., Nishimori, M., Lim, W. H., Miyazaki, C., Okamoto, A., Kanae, S., and Oki, T.: Visualizing the Interconnections Among Climate Risks, Earth's Future, 7, 85–100, 2019. 

Zaherpour, J., S. N. Gosling, N. Mount, H. M. Schmied, T. I. E. Veldkamp, R. Dankers, S. Eisner, D. Gerten, L. Gudmundsson, and I. Haddeland: Worldwide evaluation of mean and extreme runoff from six global-scale hydrological models that account for human impacts. Environ. Res. Lett., 2018. 

Zhao, M., Heinsch, F. A., Nemani, R., and Running, S. W.: Improvements of the MODIS terrestrial gross and net primary production global data set, Remote Sens. Environ., 95, 164–176,, 2005. 

Zhao, F., Veldkamp, T. I., Frieler, K., Schewe, J., Ostberg, S., Willner, S., Schauberger, B., Gosling, S. N., Schmied, H. M., Portmann, F. T., Leng, G., Huang, M., Liu, X., Tang, Q., Hanasaki, N., Biemans, H., Gerten, D., Satoh, Y., Pokhrel, Y., Stacke, T., Ciais, P., Chang, J., Ducharne, A., Guimberteau, M., Wada, Y., Kim, H., and Yamazaki, D.: The critical role of the routing scheme in simulating peak river discharge in global hydrological models, Environ. Res. Lett., 12, 075003,, 2017. 

Short summary
The most significant feature of MIROC-INTEG-LAND is that the land surface model that describes the processes of the energy and water balances, human water management, and crop growth incorporates a land-use decision-making model based on economic activities. The future simulations indicate that changes in climate have significant impacts on crop yields, land use, and irrigation water demand.