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

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 1 Relationships 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 CO₂ 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 km³ are explicitly simulated; medium-sized reservoirs with a storage capacity ranging from 3×10⁶ to 1×10⁹ m³ (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 CO₂ 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 CO₂ concentration, which is called the CO₂ fertilization effect. However, it is also known that the CO₂ 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 CO₂ fertilization effect is described as a function of atmospheric CO₂ 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 CO₂ 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 CO₂ and CH₄ 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 ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$. 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 CO₂, 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 C₃ crops such as wheat, and C₄ 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 ${\mathrm{0.5}}^{\circ }\times {\mathrm{0.5}}^{\circ }$. 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 2 The 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.

4.2 HiGWMAT + PRYSBI2

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.

4.3 VISIT

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⁻² 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; https://www.isimip.org/protocol/#isimip3b, 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 3 Comparison 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

6.1 HiGWMAT

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 4 Comparison 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 $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ 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 km³, which is greater than the 1238±67 km³ 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.

6.2 PRYSBI2

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 5 Comparison 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 6 Comparison 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).

6.3 VISIT

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⁻¹ by VISIT, 114.0 Pg C yr⁻¹ by flux upscaling, and 121.7 Pg C yr⁻¹ 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 CO₂ 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 C_j, 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 7 Comparison 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 8 Same 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 9 Same 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 10 Time 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 (km³ yr⁻¹). 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 11 Time series of changes in crop yield (unit: t ha⁻¹) 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⁻¹) 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 CO₂ fertilization effects due to rising greenhouse gas concentrations (Sect. A6) and the increase in technical coefficients due to future technological improvement (Sect. A7).

Figure 12 Time 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 CO₂ 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 13 Spatial 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 14 Temporal 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 15 Spatial 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 CO₂ concentration fixed at 2006 (noCL + noFE), those with fixed climatic factors (noCL + FE), and those with variable climatic factors and CO₂ 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 CO₂ 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 CO₂ fertilizer effects in crop models (Sakurai et al., 2014), the increase in crop yields due to the CO₂ 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 CO₂ fertilization effects (Fig. 16a). In summary, the changes in climate and CO₂ 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 CO₂ 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 16 Time series of changes in (a) cropland yield (maximum across five crops in each grid, t ha⁻¹), (b) food cropland area (a fraction of total land area), and irrigation demand (km³ yr⁻¹) based on the forcings of the five climate models under the RCP8.5 scenario. Simulations with climatic factors and CO₂ concentrations fixed at 2006 (light green, noCL + noFE), those with climatic factors fixed (cyan, noCL + FE), and those with varying climate and CO₂ 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 CO₂ 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 (l_emerge). 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 (thu_total), 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:

$$\begin{array}{}\text{(A1)}& \text{Yield}={\text{hi}}_{\mathrm{base}}\cdot {\text{BIO}}_{\mathrm{above}\left(\mathrm{maturity}\right)},\end{array}$$

where “Yield” is the crop yield (kg ha⁻¹), hi_base is the harvest index, and BIO_{above(maturity)} (kg ha⁻¹) is the aboveground biomass at the time of crop maturity. Although the harvest index changes according to atmospheric CO₂ 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 (t_d) according to the method described by Goudriaan and van Laar (1994). The amount of photosynthetically active radiation, PAR${}_{t,t_{\text{d}}}$ (MJ m⁻² s⁻¹), intercepted by the leaf at time t on a given day t_d, 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):

$$\begin{array}{}\text{(A2)}& C_{v_{c_{\text{max}}}(t,t_{\text{d}})}=\exp \left[\left({\text{TM}}_{t,t_{\text{d}}}-\mathrm{25}\right)\cdot {\displaystyle \frac{{\text{ep}}_{v_{c_{\text{max}}}}}{\mathrm{298}\cdot R\cdot (\mathrm{273}+{\text{TM}}_{t,t_{\text{d}}})}}\right],\end{array}$$

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{C}_{j_{\text{max}}(t,t_{\text{d}})}& =\exp \left[{\displaystyle \frac{E_{j_{\text{max}}}\left({\text{TM}}_{t,t_{\text{d}}}-\mathrm{25}\right)}{\mathrm{298}\cdot R\cdot \left({\text{TM}}_{t,t_{\text{d}}}+\mathrm{273}\right)}}\right]\\ & \cdot {\displaystyle \frac{\mathrm{1}+\exp \frac{\mathrm{298}\cdot S_{j_{\text{max}}}-H_{j_{\text{max}}}}{\mathrm{298}\cdot R}}{\mathrm{1}+\exp \frac{\left({\text{TM}}_{t,t_{\text{d}}}+\mathrm{273}\right)\cdot S_{j_{\text{max}}}-H_{j_{\text{max}}}}{\left({\text{TM}}_{t,t_{\text{d}}}+\mathrm{273}\right)\cdot R}}},\end{array}\end{array}$$

where $C_{v_{c_{\text{max}}}(t,t_{\text{d}})}$ and $C_{j_{\text{max}}(t,t_{\text{d}})}$ represent the effect of temperature on the maximum carboxylation rate of Rubisco and the potential rate of electron transport, respectively; TM${}_{t,t_{\text{d}}}$ is the air temperature (^∘C) at time t on day t_d; ep${}_{v_{c_{\text{max}}}}$, $E_{j_{\text{max}}}$, $S_{j_{\text{max}}}$, and $H_{j_{\text{max}}}$ 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⁻¹ K⁻¹).

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

$$\begin{array}{}\text{(A4)}& {\displaystyle }{C}_{\text{dark}(t,t_{\text{d}})}=\exp \left[\left({\text{TM}}_{t,t_{\text{d}}}-\mathrm{25}\right)\cdot {\displaystyle \frac{{\text{ep}}_{\mathrm{rd}}}{\mathrm{298}\cdot R\cdot (\mathrm{273}+{\text{TM}}_{t,t_{\text{d}}})}}\right],\end{array}$$

where C_dark(t,t_d) represents the effect of temperature on dark respiration at time t on day t_d, and ep_rd 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:

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle }{V}_{c_{\text{max}}(t,t_{\text{d}})}=\mathrm{\Theta }\cdot {\mathit{\xi }}_V\cdot C_{v_{c_{\text{max}}}(t,t_{\text{d}})}\cdot v_{c_{\text{max}}}\cdot W_{\mathrm{stress}\left(t_{\text{d}}\right)},\text{(A6)}& {\displaystyle }& {\displaystyle }{J}_{\max (t,t_{\text{d}})}=\mathrm{\Theta }\cdot {\mathit{\xi }}_J\cdot C_{j_{\text{max}}(t,t_{\text{d}})}\cdot j_{\max }\cdot W_{\mathrm{stress}\left(t_{\text{d}}\right)},\end{array}$$

where $V_{c_{\text{max}}(t,t_{\text{d}})}$ is the maximum carboxylation rate of Rubisco, $J_{max(t,t_{\text{d}})}$ is the potential rate of electron transport, and $v_{c_{\text{max}}}$ and j_max are the potential maximum carboxylation rate and the potential rate of electron transport, respectively. $W_{\mathrm{stress}\left(t_{\text{d}}\right)}$ 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 CO₂ concentration. These variables (Θ, ξ_V, and ξ_J) are described in the following section. The dark respiration rate is calculated as follows:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{R}_{\mathrm{d}(t,t_{\text{d}})}=\text{rd}\cdot C_{\mathrm{dark}\left(t,t_{\text{d}}\right)}\cdot v_{\mathrm{cmax}},\end{array}$$

where $R_{\text{d}(t,t_{\text{d}})}$ is the dark respiration rate (µmol m⁻² s⁻¹), 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 CO₂ concentration

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

$$\begin{array}{}\text{(A8)}& {\displaystyle }\begin{array}{rl}& {\mathit{\xi }}_V=\\ & \mathrm{1}-{\displaystyle \frac{r_{\mathit{\varphi }\mathrm{1}}\left(c_{\mathrm{a}}-c_{\mathrm{base}}\right)+r_{\max \mathrm{1}}-\sqrt{{\left(r_{\mathit{\varphi }\mathrm{1}}\left(c_{\mathrm{a}}-c_{\mathrm{base}}\right)+r_{\max \mathrm{1}}\right)}^{\mathrm{2}}-\mathrm{4}{r}_{\mathit{\theta }}{r}_{\mathit{\varphi }\mathrm{1}}{r}_{\text{max1}}\left(c_{\mathrm{a}}-c_{\mathrm{1}}\right)}}{\mathrm{2}{r}_{\mathit{\theta }}}},\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}& {\mathit{\xi }}_J=\\ & \mathrm{1}-{\displaystyle \frac{r_{\mathit{\varphi }\mathrm{2}}\left(c_{\mathrm{a}}-c_{\mathrm{base}}\right)+r_{\max \mathrm{2}}-\sqrt{{\left(r_{\mathit{\varphi }\mathrm{1}}\left(c_{\mathrm{a}}-c_{\mathrm{base}}\right)+r_{\max \mathrm{2}}\right)}^{\mathrm{2}}-\mathrm{4}{r}_{\mathit{\theta }}{r}_{\mathit{\varphi }\mathrm{1}}{r}_{\text{max2}}\left(c_{\mathrm{a}}-c_{\mathrm{base}}\right)}}{\mathrm{2}{r}_{\mathit{\theta }}}}\end{array},\text{(A10)}& {\displaystyle }{r}_{\mathit{\varphi }\mathrm{1}}={\displaystyle \frac{{\text{dr}}_{v_{c_{\text{max}}}}}{c_{\mathrm{base}}}},\text{(A11)}& {\displaystyle }{r}_{\mathit{\varphi }\mathrm{2}}={\displaystyle \frac{{\text{dr}}_{j_{\text{max}}}}{c_{\mathrm{base}}}},\text{(A12)}& {\displaystyle }{r}_{\max \mathrm{1}}={\text{dr}}_{v_{c_{\text{max}}}}\left({\displaystyle \frac{\mathrm{600}}{c_{\mathrm{base}}}}-\mathrm{1}\right),\text{(A13)}& {\displaystyle }{r}_{\max \mathrm{2}}={\text{dr}}_{j_{\text{max}}}\left({\displaystyle \frac{\mathrm{600}}{c_{\mathrm{base}}}}-\mathrm{1}\right),\end{array}$$

where ξ_V and ξ_J are photosynthesis compensation variables, dr${}_{v_{c_{\text{max}}}}$ and dr${}_{j_{\text{max}}}$ describe the parameters, c_a is atmospheric CO₂ concentration (mol mol⁻¹), and c_base is the baseline atmospheric CO₂ concentration (mol mol⁻¹). In this model, if dr${}_{v_{c_{\text{max}}}}$ and dr${}_{j_{\text{max}}}>\mathrm{0}$, ξ_V and ξ_J decrease linearly with increasing atmospheric CO₂. If dr${}_{v_{c_{\text{max}}}}$ and dr${}_{j_{\text{max}}}=\mathrm{0}$, ξ_V and ξ_J do not depend on atmospheric CO₂. In these equations, r_max1 and r_max2 are the respective asymptotic lines. r_θ is the parameter that determines the curvature of the lines; we set r_θ=0.99. The parameters dr${}_{v_{c_{\text{max}}}}$ and dr${}_{j_{\text{max}}}$ 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:

$$\begin{array}{}\text{(A14)}& {\displaystyle }\mathrm{\Theta }={\mathit{\theta }}_{\mathrm{base}}+{\mathit{\theta }}_{\mathrm{trend}}\left(\text{Year}-y_{\mathrm{base}}\right),\end{array}$$

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 y_base 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 y_base 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 $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ (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 ${\mathrm{30}}^{\prime \prime }\times {\mathrm{30}}^{\prime \prime }$ 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 r_i of a given 30^′′ cell being used as cropland is calculated as

$$\begin{array}{}\text{(B1)}& {\displaystyle }{r}_i={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp \left(\mathrm{1.228}+\mathrm{0.237}{\mathit{\varphi }}_i-\mathrm{0.206}{p}_k{y}_j/w_k\right)}}{C}_j,\end{array}$$

where ϕ is the slope, y is the yield per unit area (t ha⁻¹), 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 C_j 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 p_k and wage w_k 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 y_j. 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. y_j 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 w_k 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 B1 Comparison 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 B2 Global distribution of areas protected from bioenergy production.

The 0.5^∘ cell cropland area (R_j) is calculated by averaging the cropland probability in each of the 30^′′ cells (r_i) as follows:

$$\begin{array}{}\text{(B2)}& {\displaystyle }{R}_j=\sum _i^{J_i}{\displaystyle \frac{r_i}{J_i}},\end{array}$$

where J_i is the number of i cells (3600) in each 0.5^∘ cell. The adjustment parameter C_j 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 R_j 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 p_1,k and w_1,k, respectively.

First, wages in the agricultural sector, w_1,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).

$$\begin{array}{}\text{(B3)}& {\displaystyle }{l}_{\mathrm{1},k}+l_{\mathrm{2},k}=L_k,\end{array}$$

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

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\text{GDP}}_k=w_{\mathrm{1},k}\cdot l_{\mathrm{1}}+w_{\mathrm{2},k}\cdot l_{\mathrm{2}},\end{array}$$

where the value calculated by AIM/CGE is used for GDP_k (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 w_1,k is calculated as

$$\begin{array}{}\text{(B5)}& {\displaystyle }{w}_{\mathrm{1},k}={\displaystyle \frac{{\text{GDP}}_k}{l_{\mathrm{1},k}+l_{\mathrm{2},k}}}={\displaystyle \frac{{\text{GDP}}_k}{L_k}}.\end{array}$$

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 p_1,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:

$$\begin{array}{}\text{(B6)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{1},k}=w_{\mathrm{1},k}{\displaystyle \frac{l_{\mathrm{1},k}}{x_{\mathrm{1},k}}},\text{(B7)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{2},k}=w_{\mathrm{2},k}{\displaystyle \frac{l_{\mathrm{2},k}}{x_{\mathrm{2},k}}},\end{array}$$

where x_1,k and x_2,k are the production index in the agricultural and nonagricultural sectors, respectively. Here, the production index in the agricultural sector in region k (x_1,k) can be calculated as the sum of the products of 0.5^∘ crop yield y_j and cropland area R_j using Eq. (7):

$$\begin{array}{}\text{(B8)}& {\displaystyle }{x}_{\mathrm{1},k}=\sum _j^{K_j}{y}_j{R}_j,\end{array}$$

where K_j indicates the number of 0.5^∘ cells within the country/region k (3600). As described above, the cropland distribution R_j 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):

$$\begin{array}{}\text{(B9)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{1},k}={\mathit{\pi }}_k\cdot P_{\mathrm{1},k},\text{(B10)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{2},k}={\mathit{\pi }}_k\cdot P_{\mathrm{2},k}.\end{array}$$

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:

$$\begin{array}{}\text{(B11)}& {\displaystyle }{p}_{\mathrm{1},k}\cdot \left(x_{\mathrm{1},k}-X_{\mathrm{1},k}\right)+p_{\mathrm{2},k}\cdot \left(x_{\mathrm{2},k}-X_{\mathrm{2},k}\right)=\mathrm{0},\end{array}$$

where X_1,k and X_2,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 X_1,k. In this study, livestock feed demand is not included in X_1,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:

$$\begin{array}{}\text{(B12)}& {\displaystyle }{\text{GDP}}_k=P_{\mathrm{1},k}\cdot x_{\mathrm{1},k}+P_{\mathrm{2},k}\cdot x_{\mathrm{2},k}.\end{array}$$

In Eq. (B3) and Eqs. (B5)–(B11) above, the eight unknown values are p_1,k, p_2,k, x_1,k, x_2,k, l_1,k, l_2,k, π_k, and X_2,k. Of these, because the reference for the international price P is the United States (region index k=1), P_1,1 and P_2,1 (along with p_1,1, p_2,1) cannot be set. For this reason, the condition is imposed that total global net exports and imports are equal to zero:

$$\begin{array}{}\text{(B13)}& {\displaystyle }& {\displaystyle }\sum _{k=\mathrm{1}}^{K_{\text{all}}}\left(x_{\mathrm{1},k}-X_{\mathrm{1},k}\right)=\mathrm{0},\text{(B14)}& {\displaystyle }& {\displaystyle }\sum _{k=\mathrm{1}}^{K_{\text{all}}}\left(x_{\mathrm{2},k}-X_{\mathrm{2},k}\right)=\mathrm{0}.\end{array}$$

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 p_1,k, p_2,k, x_1,k, x_2,k, l_1,k, l_2,k, π_k, and X_2,k are calculated by simultaneously solving eight equations, Eq. (B3) and Eqs. (B5)–(B11), for all 17 regions ($k=\mathrm{1}-\mathrm{17}$) subject to the conditions imposed by Eqs. (B12) and (B13). The p_1,k, and w_1,k values obtained from Eq. (B4) are entered into Eq. (B1). Finally, the share of cropland for each 0.5^∘ cell R_j 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 C_j 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:

$$\begin{array}{}\text{(B15)}& {\displaystyle }{r}_{\text{bio},i}={\displaystyle \frac{C_{\text{bio}}}{\mathrm{1}+\exp \left(\mathrm{1.228}+\mathrm{0.237}{\mathit{\varphi }}_i-\mathrm{0.206}{p}_{\text{bio},k}{y}_{\text{bio},j}/w_{\mathrm{1},k}\right)}},\end{array}$$

where ϕ_i is the slope in 30^′′ cell i, p_bio,k is the biofuel crop price in region k, y_bio,j is the yield (t ha⁻¹) of biofuel crops in 0.5^∘ cells, and w_1,k is the agricultural sector wage in region k. For the biofuel crop price p_bio,k, the values generated by AIM/CGE are used. For biofuel crop yield y_bio,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⁻¹ yr⁻¹ if the fertilizer input at the country level was below 160 kg(N) ha⁻¹ yr⁻¹ 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 C_bio 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:

$$\begin{array}{}\text{(B16)}& \sum _k^{K_{\text{all}}}{X}_{\text{bio},k}=\sum _j^{J_{\text{all}}}{y}_{\text{bio},j}{R}_{\text{bio},j},\end{array}$$

where X_bio,k is the biofuel crop demand for region k calculated by AIM/CGE and K_all and J_all are the total number of regions (17) and the total number of 0.5^∘ cells (259 200), respectively. R_bio,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×10⁶ km², and the KBA covered an area of 19.9×10⁶ km². 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

$$\begin{array}{}\text{(B17)}& r_{\text{past},i}={\displaystyle \frac{C_{\text{past},j}\times {\text{NPP}}_j}{\left(\mathrm{1}+\mathit{\varphi }/\mathrm{20}\right)}}.\end{array}$$

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 NPP_j. The boundary condition of the VISIT offline simulations is fixed at year 2005. C_past,j is the adjustment parameter for 0.5^∘ cells. The value of C_past,j changes from year to year. The adjustment parameter for the base year, $C_{\text{past},j}(t=\mathrm{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, C_past,j(t), are set by applying a proportionality factor, α(t), to the base-year parameter:

$$\begin{array}{}\text{(B18)}& C_{\text{past},j}\left(t\right)=\mathit{\alpha }\left(t\right)\times C_{\text{past},j}(t=\mathrm{0}),\end{array}$$

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

$$\begin{array}{}\text{(B19)}& S_{\text{past},k}\left(t\right)=\sum _j^{J_k}{R}_{\text{past},j}\left(t\right)\end{array}$$

is met, where S_past,k(t) is the pastureland demand calculated by AIM/CGE for region k; R_past,j(t) is the average of percentage of pastureland for all 30^′′ cells (from Eq. B16) in a given 0.5^∘ cell; and J_k 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, R_mfr,j(t), is formularized in terms of the area of managed forests in the base year and the population density:

$$\begin{array}{}\text{(B20)}& R_{\text{mfr},j}=A_{\text{fr},j}\times {\displaystyle \frac{{\mathit{\rho }}_{j^{*}}}{C_{\text{manfr},k}+{\mathit{\rho }}_{j^{*}}}},\end{array}$$

where A_fr,j is the area of managed forest in the 0.5^∘ cells in the base year (2005) and ${\mathit{\rho }}_{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⁻¹ t⁻¹) is balanced with the price of timber (∼ USD 100 t⁻¹). 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 (${\mathit{\rho }}_{j^{*}}$). C_mfr,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 X_mfr,k (kg yr⁻¹) calculated by AIM/CGE. We use the region-level adjustment factors for managed forest (C_mfr,k), because the grid-level reference data are not available. In other words, C_mfr,k is set so that the total regional amount of harvested timber equals the regional total roundwood demand:

$$\begin{array}{}\text{(B21)}& X_{\text{mfr},k}=\sum _j^{J_k}{R}_{\text{mfr},j}\times {\displaystyle \frac{B_j}{L_j}},\end{array}$$

where B_j is the distribution of forest biomass (kg m⁻²) 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). J_k is the total number of 0.5^∘ cells in each region k. L_j is the harvesting period (years), which is estimated as follows, based on theNPP_j for 0.5^∘ cells obtained from VISIT (Ito and Inatomi, 2012).

$$\begin{array}{}\text{(B22)}& L_j=\left\{\begin{array}{ll}\mathrm{\infty }& {\text{NPP}}_j<\mathrm{4}\\ \mathrm{500}/{\text{NPP}}_j& \mathrm{4}\le {\text{NPP}}_j\le \mathrm{25}\\ \mathrm{20}& \mathrm{25}<{\text{NPP}}_j\end{array}\right.\end{array}$$

L_j 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 $R_{\text{mfr},j}\times B_j/L_j$ (kg yr⁻¹) based on the distribution of managed forests R_mfr,j, forest biomass B_j, and the felling period L_j 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 B1 Correspondence of land-cover type in land-use model and transition matrix.

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 B3 Schematic 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.