Assessing effects of climate and technology uncertainties in large natural resource allocation problems

The productivity of the world's natural resources is critically dependent on a variety of highly uncertain factors, which obscure individual investors and governments that seek to make long-term, sometimes irreversible, investments in their exploration and utilization. These dynamic considerations are poorly represented in disaggregated resource models, as incorporating uncertainty into large-dimensional problems presents a challenging computational task. In this paper, we apply the SCEQ algorithm (Cai and Judd, 2023) to solve a large-scale dynamic stochastic global land resource use problem with stochastic crop yields due to adverse climate impacts and limits on further technological progress. For the same model parameters and bounded shocks, the range of land conversion is considerably smaller for the dynamic stochastic model than for deterministic scenario analysis.

1 Introduction

Understanding the future allocation of the world's natural resources is an important research problem for environmental scientists and economists. This involves a thorough understanding of a complex interplay of different factors, including, among others, continuing population increases, shifting diets among the poorest populations in the world, increasing production of renewable energy, including biofuels, and growing demand for ecosystem services, including forest carbon sequestration (Foley et al., 2011). The problem is further confounded by faster-than-expected climate change which is altering the biophysical environment of agriculture and forestry. Moreover, highly uncertain future productivities and valuations of ecosystem services, coupled with medium- to long-term irreversibilities in the extraction of nonrenewable or partially renewable resources, such as long-growth natural forests, give rise to a challenging problem of decision-making under uncertainty.

While there is a large body of research analyzing the problem of natural resource extraction and utilization under uncertainty theoretically or using stylized computational models (see, e.g., Miranda and Fackler, 2004; Tsur and Zemel, 2014, and references therein), quantifying the effects of uncertainty in the natural resource use in a more realistic setting remains a challenging problem. This is because natural resource allocation problems, like environmental policy problems in general, involve highly nonlinear structure and damage functions, important irreversibilities, and long time horizons (Pindyck, 2007). Computationally integrated models of economy and environment are the standard workhorse mechanisms for modeling the long-term allocation of the world's natural resources, including particularly difficult land use problems (see, e.g., Füssel, 2009; Schmitz et al., 2014; Nikas et al., 2019, and references therein). These models have the important advantage of detailed spatial and sectoral (particularly in the energy and agricultural sector) coverage, which allows them to capture a broad range of responses to changes in demand and supply factors affecting the utilization of natural resources. However, given the high computational complexity of these models, they are typically either static or based on myopic expectations, whereby decisions about production, consumption, and resource extraction and conversion are made only on the basis of information in the period of the decision (Babiker et al., 2009). These models, therefore, have limited ability to address important intertemporal questions such as, for example, a dynamic trade-off between conservation, carbon sequestration, and renewable offsets for fossil fuels. Among the few forward-looking dynamic economy and environment models, none explicitly incorporates uncertainty into the determination of the optimal path of natural resource use.¹ This is because introducing uncertainty into these models is confined by an array of computational obstacles that are very difficult (e.g., high dimensionality and kinks caused by occasionally binding constraints), if not impossible, to address using standard numerical methods such as projection methods and value function iteration (see, e.g., Judd, 1998; Miranda and Fackler, 2004; Cai and Judd, 2014; Cai, 2019). To the extent that uncertainty in these models is considered, this is only through parametric or probabilistic sensitivity analysis or the use of alternative scenarios. Therefore, the high-dimensional resource use models have not effectively dealt with optimal extraction and conversion decisions along the uncertain path of key drivers affecting resource allocation in the face of costly reversal of conversion decisions.

In this study, we seek to address this important limitation of the economy–environment modeling of natural resource use. In doing so, we build on recent advances in computational economic and operation research. Cai et al. (2017) introduced a nonlinear certainty equivalent approximation method (NLCEQ) for solving large-scale, infinite-horizon stationary dynamic stochastic problems and demonstrated how this method could be used to achieve an accurate solution to a stylized stationary dynamic stochastic land use problem. The original NLCEQ method, however, is ill-suited for solving most environmental and resource economic problems. This is because stochastic problems of the utilization of natural resources feature nonstationary stochastic trends, such as climate or technological trajectories, and some never converge to a stationary state. Building on the original NLCEQ work Cai and Judd (2023) introduced a simulated certainty equivalent approximation method (SCEQ), which efficiently solves nonstationary dynamic stochastic problems, including those with high dimensionality and occasionally binding constraints. Cai and Judd (2023) showed that the SCEQ method is highly accurate and achieves stable numerical solutions for dynamic stochastic problems in economics.

We apply the SCEQ method to solve a large-scale dynamic stochastic model, focusing on the optimal global land use allocation problem. This highly complicated resource use problem features multiple cross-sectoral and dynamic trade-offs. Specifically, we apply the method to a global land use model nicknamed FABLE (Forest, Agriculture, and Biofuels in a Land use model with Environmental services) in the face of uncertainty. FABLE is a dynamic, forward-looking, global, multi-sectoral, and partial-equilibrium model designed to analyze the evolution of global land use over the coming century. Prior applications of that model (Steinbuks and Hertel, 2013; Hertel et al., 2013, 2016; Steinbuks and Hertel, 2016) analyze the competition for scarce global land resources in light of the growing demand for food, energy, forestry, and environmental services and evaluate key drivers and policies affecting global land use allocation. All of these applications, however, assume perfect foresight and treat uncertainty in a parametric fashion, thus ignoring the impact of future uncertainties in the optimal allocation of global land use. To ensure consistency between theoretical and numerical model solutions, we assume the bounded solution space. As we show below, this assumption is well justified for economic models of large natural resource allocation problems, including the FABLE model.

By way of illustration, we focus on uncertainty emanating from crop productivity over the next century. Along with energy prices, regulatory policies, and technological change in food, timber, and biofuels industries, this is one of four core uncertainties affecting competition for global land use (Steinbuks and Hertel, 2013). To quantify the uncertainty in agricultural yields, we construct a stochastic crop productivity index that captures two key uncertainty sources: technological progress and global climate change (Lobell et al., 2009; Licker et al., 2010; Foley et al., 2011).² Following Rosenzweig et al. (2014), we use projections from climate and crop-simulation models under Representative Concentration Pathways 6.0 W m⁻² (RCP6) greenhouse gas (GHG)-forcing scenario (Moss et al., 2008), as well as the survey of recent agro-economic and biophysical studies to calibrate the index.

We simulate the results of the model, where the global planner optimally allocates land uses under the perfect foresight of different realizations of the crop productivity index, focusing our attention on the current century. We then compare and contrast them with the results of the dynamic stochastic model, where the global planner has rational expectations about crop yields subject to bounded autocorrelated climate shocks. When the uncertainty in crop productivity is incorporated into the model, we see an additional redistribution of land resources to offset the impact of potentially lower yields. Due to intertemporal substitution, some of that redistribution occurs even in the absence of actual changes in the states of climate or technology affecting crop yields. Moreover, the range of these alternative optimal paths of cropland is considerably smaller than the magnitude of possible land conversion resulting from the scenario analysis based on the deterministic model. This result indicates that when the climate shocks are bounded, the scenario analysis may significantly overstate the expected agricultural land conversion magnitude under uncertain crop yields.

Our study contributes to the growing literature that analyzes the intertemporal allocation of land and other natural resources under uncertainty and irreversibility constraints. Most of the literature focuses on a particular type of resource or sector in which intertemporal issues are significant and cannot be ignored. One example of this literature is forestry management in the context of uncertain fire risks and climate mitigation policies (Sohngen and Mendelsohn, 2003, 2007; Daigneault et al., 2010). Another example is natural land conservation decisions under irreversible biodiversity losses (Conrad, 1997, 2000; Bulte et al., 2002; Leroux et al., 2009). While these models are undoubtedly helpful for understanding the broad implications of uncertainty in the intertemporal allocation of land resources, they are less effective in quantifying the effect of uncertainty in the supply and demand drivers in more complex settings, such as, e.g., the optimal allocation of multiple competing land resources in the long run.

Our study is perhaps most closely related to the recent works of Lanz et al. (2017) and Zhao et al. (2021). Lanz et al. (2017) developed a two-sector stochastic Schumpeterian growth model with the endogenous allocation of global land use. As our paper does, they find that the optimal allocation of global land use requires more cropland conversion when the uncertainty in agricultural productivity is present. Lanz et al. (2017), however, focused on endogenous population dynamics, labor allocation, and technological progress, whereas our paper is concerned with the endogenous allocation of multiple types of land use and corresponding land-based goods and services. Our paper also advances the methodological grounds by applying a more advanced algorithm that overcomes computational difficulties in solving multidimensional stochastic land use models, which made Lanz et al. (2017) significantly simplify their model by assuming that their binary shocks occur only in three time periods. Zhao et al. (2021) compare models with adaptative expectations and perfect foresight assumptions for both the price and yield for agricultural producers to make land allocation and production decisions. Zhao et al. (2021) find similar results that land use change variation becomes much smaller than in the perfect foresight model, which allows for faster land use adjustments while market price variations increase. Unlike our paper, Zhao et al. (2021) do not explicitly incorporate uncertainty in the model's optimization stage.

2 Stochastic FABLE model

This section presents a modeling framework for analyzing nonlinear dynamic stochastic models of natural resource use with multiple sectors, in which preferences, production technology, resource endowments, and other exogenous state variables evolve stochastically over time according to a Markov process with time-varying transition probabilities. The constructed model belongs to the class of stochastic growth models, with multiple sectors studied in Brock and Majumdar (1978), Majumdar and Radner (1983), and Stokey et al. (1989) among others. Similar to other models in this class, the FABLE model has a bounded solution space, as all these models are theoretically shown to have equilibrium paths (Stokey et al., 1989). This assumption is important because it is often impossible to prove that in the presence of unbounded solution space, the stochastic model has a finite solution. For example, Weitzman's dismal theorem (Weitzman, 2010) shows that a fat-tail damage function with an infinite upper bound leads to an infinite risk premium, but a numerical truncation to finite support will always have a finite risk premium. So, assuming the bounded solution space is necessary for avoiding the potential qualitative inconsistencies between their theoretical and numerical results. We further discuss the assumption of the bounded solution space in the FABLE model in Sect. 5.

Specifically, we develop a stochastic version of a global land use model nicknamed FABLE (Forest, Agriculture, and Biofuels in a Land use model with Environmental services), a dynamic multi-sectoral model for the world's land resources over the next century (Steinbuks and Hertel, 2012, 2016). This model combines recent strands of agronomic, economic, and biophysical literature into a single intertemporally consistent analytical framework at the global scale. FABLE is a discrete dynamic partial-equilibrium model where the population, labor, physical and human capital, and other variable inputs are assumed to be exogenous. Total factor productivity and technological progress in non-land-intensive sectors are also predetermined. The model focuses on the optimal allocation of scarce land across competing uses across time and solves the dynamic paths of alternative land uses, which together maximize global economic welfare.

Figure 1 Structure of the FABLE Model. Note that state variables are shown as oval shapes. Decision variables are shown as rectangular shapes. The utility function is shown as an octagonal shape. Stochastic model terms incorporating random processes are shown as dashed shapes or arrows.

The FABLE model accommodates a complex dynamic interplay between different types of global land use, whereby the societal objective function places value on processed crops and livestock, energy services, timber products, ecosystem services, and other non-land goods and services (Fig. 1). There are three accessible primary resources in this partial-equilibrium model of the global economy: land, liquid fossil fuels, and other primary inputs, e.g., labor and capital (see the bottom part of Fig. 1). The supply of land is fixed and faces competing uses that are determined endogenously by the model. They include unmanaged forest lands which are in an undisturbed state (e.g., parts of the Amazon tropical rainforest ecosystem), agricultural (or crop) land, pasture land, and commercially managed forest land. As trees of different ages have different timber yields and different propensities to sequester carbon, the model keeps track of various tree vintages in managed forests, which introduces additional numerical complexity for solving the model. We do not keep track of vintages for natural lands and assume they are primarily old-grown forests. We ignore other land use types, such as savannah, grasslands, and shrublands, which are largely unmanaged and often of limited productivity.³ We also ignore residential, retail, and industrial uses of land in this partial-equilibrium model of agriculture and forestry.

The flow of liquid fossil fuels evolves endogenously along an optimal extraction path, allowing exogenously specified new discoveries of fossil fuel reserves. Other primary inputs include variable inputs, such as labor, capital (both physical and human), and intermediate materials. The endowment of other primary inputs is exogenous and evolves along a pre-specified global economic growth path.

There are six intermediate inputs used in the production of land-based goods and services in FABLE: petroleum products, fertilizers, crops, liquid biofuels,⁴ live animals, and raw timber (see the middle part of Fig. 1). Fossil fuels are refined and converted to either petroleum products that are further combusted or to fertilizers that are used to boost yields in the agricultural sector. Cropland and fertilizers are combined to grow crops that can be further converted into processed food and biofuels or used as animal feed. Specifically, we assume that agricultural land $L_t^{\mathrm{A},\mathrm{c}}$ and fertilizers $x_t^{n,\mathrm{c}}$ are imperfect substitutes in the production of food crops, $x_t^{\mathrm{c}}$, with specific production technology given by the following constant elasticity of substitution (CES) function:

$$\begin{array}{}\text{(1)}& x_t^{\mathrm{c}}={\mathit{\theta }}_t^{\mathrm{c}}{\left({\mathit{\alpha }}^n{\left(L_t^{\mathrm{A},\mathrm{c}}\right)}^{{\mathit{\rho }}_n}+\left(\mathrm{1}-{\mathit{\alpha }}^n\right){\left(x_t^{n,\mathrm{c}}\right)}^{{\mathit{\rho }}_n}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_n}},\end{array}$$

where ${\mathit{\theta }}_t^{\mathrm{c}}$ is stochastic crop technology index, and αⁿ and ρ_n are, respectively, the input share and substitution parameters. Equation (1) captures three key responses within the model to changes in crop technology index: (i) demand response (change in consumption of food crops), (ii) adaptation on the extensive margin (substitution of agricultural land for other land resources), and (iii) adaptation on the intensive margin (substitution of agricultural land for fertilizers).

Biofuels substitute imperfectly for liquid fossil fuels in final energy demand. The food crops used as animal feed and pasture land are combined to produce raw livestock. Harvesting managed forests yields raw timber that is further used in timber processing.

The land-based consumption goods and services take the form of processed crops, livestock, and timber and are, respectively, outcomes of food crops, raw livestock, and timber processing. The production of energy services combines non-land energy inputs (i.e., liquid fossil fuels) with biofuels, and the resulting mix is further combusted. Finally, all land types can contribute to other ecosystem services and a public good to society, including recreation, biodiversity, and other environmental goods and services. To close the demand system, we also include other non-land goods and services (e.g., manufacturing goods and retail, construction, financial, and information services), which involve the “consumption” of other primary inputs not spent on the production of land-based goods and services. As the model focuses on the representative agent's behavior, the final consumption products are all expressed in per capita terms.

A complete description of model equations, variables, and parameter values is presented in the Appendix.

3 Modeling crop yield uncertainty

This section characterizes uncertainty in future agricultural yields over the coming century, which is one of the core uncertainties shown to affect land use in the long run (Steinbuks and Hertel, 2013). Crop yields are subject to two types of uncertainties: those related to the development and dissemination of new technologies and those related to changes in the climatic conditions under which the crops are grown. The former type of uncertainty has until recently dominated the pattern of the evolution of global crop yields, whereas the latter is becoming an increasingly important factor (Lobell and Field, 2007; IPCC, 2014). While it is plausible to hypothesize that accelerating climate impacts may, in turn, induce further technological advances in an effort to facilitate adaptation to climate change, this hypothesis is not supported by limited empirical evidence (Burke and Emerick, 2016). Therefore, in this paper, these two sources of uncertainty are treated separately, although they are both characterized by the use of combined climate and crop-simulation models run over a global grid.

We characterize future uncertainty in yields by constructing a stochastic crop productivity index, ${\mathit{\theta }}_t^{\mathrm{c}}$, which captures the evolution of future crop yields under different realizations of uncertainty in crop productivity based on the most recent projections in the agronomic and environmental science studies. An important characteristic of staple grain yields is that they tend to grow linearly, adding a constant amount of gain (e.g., t ha⁻¹) each year (Grassini et al., 2013). This suggests that the proportional growth rate should fall gradually over time. However, crop physiology dictates certain biophysical limits to the rate at which sunlight and soil nutrients can be converted to grain. And there is some recent agronomic evidence (Cassman et al., 2010; Grassini et al., 2013) showing that yields appear to be reaching a plateau in some of the world's most important cereal-producing countries. Cassman (1999) suggests that average national yields can be expected to plateau when they reach 70 %–80 % of the genetic yield potential ceiling. Based on these observations from the agronomic literature, we specify the following logistic function determining the evolution of the crop productivity index over time:

$$\begin{array}{}\text{(2)}& {\mathit{\theta }}_t^{\mathrm{c}}={\displaystyle \frac{{\mathit{\theta }}_T^{\mathrm{c}}{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{c}}{e}^{{\mathit{\kappa }}_{\mathrm{c}}t}}{{\mathit{\theta }}_T^{\mathrm{c}}+{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{c}}\left(e^{{\mathit{\kappa }}_{\mathrm{c}}t}-\mathrm{1}\right)}},\end{array}$$

where ${\mathit{\theta }}_{\mathrm{0}}^{\mathrm{c}}$ is the value of the crop productivity index in period 0, which we calibrate to match observed weighted yields in key staple crops (corn, rice, soybeans, and wheat). ${\mathit{\theta }}_T^{\mathrm{c}}$ is the crop yield potential at the end of the current century; that is, “the yield an adapted crop cultivar can achieve when crop management alleviates all abiotic and biotic stresses through optimal crop and soil management” (Evans and Fischer, 1999). κ_c is the logistic convergence rate to achieving potential crop yields.

Though the initial value of the crop productivity index is known with certainty, potential crop yields are highly uncertain. We assume that potential crop yields are affected by a two-dimensional stochastic process of climate and technological shocks, J_1,t, and J_2,t, respectively. For the technological shock, J_2,t, we assume that there are three states of technology: “bad” (indexed by $J_{\mathrm{2},t}=\mathrm{1}$), “medium” (indexed by $J_{\mathrm{2},t}=\mathrm{2}$), and “good” (indexed by $J_{\mathrm{2},t}=\mathrm{3}$). In the optimistic (i.e., good) state of advances in crop technology, we assume that yields will continue to grow linearly throughout the coming century, eliminating the yield gap by 2100. In the medium state of technology, rather than closing the yield gap by 2100, average yields in 2100 are just three-quarters of the yield potential at that point in time. In the bad state of technology, there is no technological progress, and the crop yields stay the same as at the beginning of the coming century.

For the climate shock, J_1,t, we assume it is a Markov chain with five possible states at each time t. To construct these states, we use the results of Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) fast-track crop-simulation model comparison (Rosenzweig et al., 2014).⁵ Since FABLE is a partial-equilibrium model without an embedded climate module, we cannot directly capture all sources of GHG emissions and endogenize their effect on global crop yields. Instead, to ensure the simulation results' comparability with the structural parameters (e.g., demographic and economic growth and the rate of technological change) of the FABLE model, we select four crop-simulation model runs under the RCP6 GHG forcing scenario.⁶ We also consider alternative assumptions on CO₂ fertilization effects.

Based on these data, we construct five states corresponding to quintiles of the distribution of different outcomes of four global crop-simulation models and five global climate models, with and without CO₂ fertilization effects for potential crop yields by 2100. Under two optimistic states of the world, we observe 2 % and 15 % increases in potential crop yields relative to the model baseline (calibrated based on historical trends over the reference period 1971 to 2004), respectively, whereby significant CO₂ fertilization effects offset the negative effects of climate change. For the next two states, we see declines of 15 % and 19 % in potential crop yields relative to the model baseline, whereby CO₂ fertilization effects are assumed to be either small or non-existent, and the negative effects of climate change tend to prevail. Finally, under the most pessimistic state of the world, drastic adverse effects of climate change combined with the absence of any CO₂ fertilization effects result in a 36 % decline in potential crop yields relative to the model baseline.

Further details of constructing climate and technological states can be found in the Appendix.

The path of technological change in crop yields evolves by reversible transitions across these states. The stochastic path of the crop productivity index is then given by

$$\begin{array}{}\text{(3)}& A_t={\displaystyle \frac{A_T(J_{\mathrm{1},t},J_{\mathrm{2},t})A_{\mathrm{0}}{e}^{{\mathit{\kappa }}_{\mathrm{c}}t}}{A_T(J_{\mathrm{1},t},J_{\mathrm{2},t})+A_{\mathrm{0}}\left(e^{{\mathit{\kappa }}_{\mathrm{c}}t}-\mathrm{1}\right)}},\end{array}$$

where A_T(i,j) represents the crop productivity index at the terminal time T at the state $J_{\mathrm{1},t}=i$ and $J_{\mathrm{2},t}=j$ for $i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{5}$ and $j=\mathrm{1},\mathrm{2},\mathrm{3}$. Thus, A_t is a Markov chain, which takes 1 of 15 possible time-varying values at each time period. This can be seen as a discretization of a mean-reverting process with continuous values and a time trend, but a finer Markov chain with more values can only marginally change our solution. As A_t is completely dependent on J_1,t and J_2,t, it is not a state variable, whereas J_1,t and J_2,t are both state variables.

Having characterized the realizations of crop productivity under alternative states of agricultural technology and climate change impacts, we still need to calibrate the transition probabilities for the climate and technology shocks to construct the stochastic crop productivity index. As regards climate shock, the environmental and climate science literature acknowledges some degree of persistence but does not provide much guidance on the transition dynamics between alternative climate states affecting crop yields. In the absence of reliable estimates for constructing the transition probability matrix of J_1,t, we assume simple transition dynamics, where each state has a 50 % probability of retaining itself next period and a 25 % probability of moving upwards and downwards to an adjacent state (note that realizations can only stay the same or move upwards from the lowest state, e.g., $J_{\mathrm{1},t}=\mathrm{1}$, and only stay the same or move downwards from the highest state, e.g., $J_{\mathrm{1},t}=\mathrm{5}$). As regards the technology shock, since we do not have historical data on the evolution of agricultural technology, we assume that technological advances in agriculture follow a similar trend to advances in the rest of the economy and use the probability transition matrix of J_2,t estimated by Tsionas and Kumbhakar (2004) for a comprehensive panel of 59 countries over the period of 1965–1990. These estimates correspond to a 20 % probability of the bad technological state, 56 % of the medium state, and 24 % of the good state. The transition probability matrices of J_1,t and J_2,t are shown in the Appendix.

Figure 2 Crop productivity index.

Figure 2 shows the deterministic-baseline path (the solid line) used in the perfect foresight model and the range of the stochastic crop productivity index based on 1000 simulation paths over the entire 21st century, with additional summary statistics presented in the Appendix. The simulations start at the medium states of climate and technology in the initial year. The deterministic-baseline path is calculated by taking expectations of the stochastic crop productivity index conditional on the initial medium states (Eq. 3). It also takes the same values as the median path (the “o” line) of simulations, whereby the climate and technological states are kept at medium, while the average line (the “+” line) deviates a bit after the year 2070. At every time t, there are 1000 realized values of A_t, among which there are only 15 different values. The 10 % and 90 % quantile lines (the dashed and dash-dotted lines) represent the 10 % and 90 % quantiles of these 1000 simulated values of A_t at time t, so they are not the realized sample paths, but Fig. 2 also displays one realized sample path of A_t, which is the dotted line.

4 Method of model solution: the SCEQ algorithm

In most resource use problems under uncertainty, including the FABLE model, the social planner's problem cannot be solved analytically, although certain inferences about the potential effects of uncertainty can be made from more stylized models. Numerical dynamic programming with value function iteration (see, e.g., Cai and Judd, 2014; Cai, 2019) is a typical method to solve these dynamic stochastic problems. However, numerical dynamic programming faces challenging problems such as high dimensionality of state space, shape-preservation of value functions (Cai and Judd, 2013), and kinks caused by occasionally binding constraints. These challenges are common in modeling natural resource use and are hard to address, even with the most advanced methods such as parallel dynamic programming (Cai et al., 2015).

For non-stationary problems, value function iteration involves computing decision rules at each period t. However, computing all these rules can be very time-consuming and unnecessary if our primary goal is to obtain simulation paths and their distributions until a time of interest, T^*, even if it is a long one (in environmental and climate change economics, for example, we are often interested in solutions for the coming century and set the time of interest to 100 years and the problem horizon of more than 300 years to avoid a large impact of terminal conditions). Instead of solving for optimal decisions for all possible states at each time, we can approximately solve for optimal decisions for those simulated states along simulated paths. This logic is embedded in the novel SCEQ algorithm (Cai and Judd, 2023) used to solve for a simulated range of land use trajectories in the stochastic FABLE model.

Below we present the SCEQ algorithm version for solving the finite horizon nonstationary stochastic dynamic programming problems that the FABLE model belongs to (for a detailed description of other cases, see Cai and Judd, 2023). Following the standard notation in the literature, let S_t be a vector of state variables (here land cover types and stocks of fossil fuels), and let a_t be a vector of decision variables (here land conversion, resource extraction, transformation, and final consumption of land-based goods and services) at each time t. The transition law of the state vector S is

$${\mathit{S}}_{t+\mathrm{1}}={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t,{\mathit{ϵ}}_t),$$

where ϵ_t is a serially uncorrelated random vector process, and G_t is a vector of functions; its ith element, G_t,i, returns the ith state variable at t+1: $S_{t+\mathrm{1},i}$. For simplicity, we assume the mean of ϵ_t is 0.⁷

The FABLE model equations (see Appendix) can be compactly represented by the following social planner's problem:

$$\begin{array}{}\text{(4)}& \begin{array}{rlr}\underset{{\mathit{a}}_t}{max}& & \mathbb{E}\left\{\sum _{t=\mathrm{0}}^{T-\mathrm{1}}{\mathit{\delta }}^t{\mathcal{U}}_t\left({\mathit{S}}_t,{\mathit{a}}_t\right)+{\mathit{\delta }}^T{V}_T\left({\mathit{S}}_T\right)\right\}\\ \mathrm{s}.\mathrm{t}.& & {\mathit{S}}_{t+\mathrm{1}}={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t,{\mathit{ϵ}}_t),t=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },T-\mathrm{1},\\ & & {\mathbf{F}}_t({\mathit{S}}_t,{\mathit{a}}_t)\ge \mathrm{0},t=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },T-\mathrm{1},\end{array}\end{array}$$

where 𝒰_t is a utility function, $\mathit{\delta }\in \left(\mathrm{0},\mathrm{1}\right)$ is the discount factor, 𝔼 is the expectation operator, T is the horizon (T=∞ if it is an infinite-horizon problem), V_T(S_T) is a given terminal value function depending on the terminal state S_T (it is zero everywhere for an infinite-horizon problem), and ${\mathbf{F}}_t({\mathit{S}}_t,{\mathit{a}}_t)\ge \mathrm{0}$ is a vector of feasibility constraints of actions a_t at time t. And we assume that the initial state S₀ is given, as it can usually be observed or estimated.

Algorithm 1 SCEQ for finite-horizon stochastic dynamic programming problems with time-variant exogenous paths.

Step 1.

Initialization step. Given the initial state S₀ and a time of interest T^*, as well as a terminal value function V_T(S_T), simulate a sequence of ϵ_t to get m paths, denoted ${\mathit{ϵ}}_t^i$ for path i, from t=0 to $T^{*}-\mathrm{1}$. Let ${\mathit{S}}_{\mathrm{0}}^i={\mathit{S}}_{\mathrm{0}}$ and iterate forward through steps 2 and 3 for $s=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },T^{*}-\mathrm{1}$.

Step 2.

Optimization step. Solve the following deterministic model starting from time s and simulated node ${\mathit{S}}_s^i$:

$$\begin{array}{}\text{(5)}& \begin{array}{rlr}\underset{{\mathit{a}}_t}{max}& & \sum _{t=s}^{T-\mathrm{1}}{\mathit{\delta }}^{t-s}{\mathcal{U}}_t\left({\mathit{S}}_t,{\mathit{a}}_t\right)+{\mathit{\delta }}^{T-s}{V}_T\left({\mathit{S}}_T\right)\\ \mathrm{s}.\mathrm{t}.& & {\mathit{S}}_{t+\mathrm{1}}={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t,\mathrm{0}),t=s,s+\mathrm{1},\mathrm{\dots },T-\mathrm{1},\\ & & {\mathbf{F}}_t({\mathit{S}}_t,{\mathit{a}}_t)\ge \mathrm{0},t=s,s+\mathrm{1},\mathrm{\dots },T-\mathrm{1},\end{array}\end{array}$$

where S_s is given by ${\mathit{S}}_s^i$, for each $i=\mathrm{1},\mathrm{\dots },m$.

Step 3.

Simulation step. Set ${\mathit{S}}_{s+\mathrm{1}}^i={\mathit{G}}_t({\mathit{S}}_s^i,{\mathit{a}}_s^i,{\mathit{ϵ}}_s^i)$, where ${\mathit{a}}_s^i$ is the optimal decision at time s of problem (5) for each $i=\mathrm{1},\mathrm{\dots },m$.

Algorithm 1 obtains simulated pathways of optimal decisions and states. Note that the inside loop across i can be switched with the outside loop across time; that is, for each i, we can obtain one simulation path by iteratively solving Eq. (5) and simulating ${\mathit{S}}_{s+\mathrm{1}}^i={\mathit{G}}_t({\mathit{S}}_s^i,{\mathit{a}}_s^i,{\mathit{ϵ}}_s^i)$ for $s=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },T^{*}-\mathrm{1}$.

The optimization step of Algorithm 1 applies the certainty equivalent approximation idea of the NLCEQ method (Cai et al., 2017); for a given state at time s, ${\mathit{S}}_s^i$, we replace all future stochastic variables by their corresponding expectations conditional on the current state ${\mathit{S}}_s^i$,⁸ and convert the dynamic stochastic problem (4) into a deterministic finite-horizon dynamic problem (5).

We implement the optimal control method (Cai, 2019) to solve Eq. (5) numerically; that is, we view Eq. (5) as a large-scale nonlinear constrained optimization problem with $\mathit{\{}{\mathit{a}}_t^i:t\ge s\mathit{\}}$ and $\mathit{\{}{\mathit{S}}_t^i:t\ge s\mathit{\}}$ as its variables and the transition equations and feasibility restrictions as its constraints. The problem can be directly solved with an appropriate nonlinear optimization solver such as CONOPT (Drud, 1994).

Observe that we just need to save the solution of Eq. (5) at time s, ${\mathit{a}}_s^i$, for use in the next step. In step 3 of Algorithm 1, we use the saved optimal decision ${\mathit{a}}_s^i$ to generate the next-period state, ${\mathit{S}}_{s+\mathrm{1}}^i={\mathit{G}}_t({\mathit{S}}_s^i,{\mathit{a}}_s^i,{\mathit{ϵ}}_s^i)$, given the realization of shocks, ${\mathit{ϵ}}_s^i$. Once we reach the state ${\mathit{S}}_{s+\mathrm{1}}^i$ at time s+1, we come back to implement step 2 and then step 3. In other words, Algorithm 1 uses an adaptive management way; decisions are made for the current period in the face of future uncertain shocks. Once the next-period shock is observed, decisions for the next period are made with re-optimization, given the observed shock and new state variables at the next period. Observe that the serial correlation of random variables has been captured in their associated transition laws. Repeating this process iteratively through T^* times, we compute a representative simulated pathway of optimal decisions, $\mathit{\{}{\mathit{a}}_s^i{\mathit{\}}}_{s=\mathrm{0}}^{T^{*}-\mathrm{1}}$, and states, $\mathit{\{}{\mathit{S}}_s^i{\mathit{\}}}_{s=\mathrm{0}}^{T^{*}}$, which corresponds to the realized path of shocks, $\mathit{\{}{\mathit{ϵ}}_s^i{\mathit{\}}}_{s=\mathrm{0}}^{T^{*}-\mathrm{1}}$. Repeating over i, we compute m simulated paths of optimal states and decisions and then obtain their distributions. This simulation process can be naturally parallelized.

After we obtain the simulated solutions for our dynamic stochastic land use problem, we also check the normalized Euler errors and find that the ℒ¹ error in the solutions for the first 100 years (the periods of interest) among 1000 simulated paths is only $\mathrm{8.6}\times {\mathrm{10}}^{-\mathrm{4}}$, and the corresponding ℒ^∞ error is only 0.02. This is within range of acceptable accuracy for the most dynamic stochastic problems (Cai and Judd, 2023).

5 Model results

Below, we describe the results of the impact of crop yield uncertainty in the optimal path of global land use based on the dynamic stochastic model simulations. We solve the model over 400 years with 5-year time steps and present the results for the first 100 years to minimize the effect of terminal period conditions on our analysis.⁹ We first present the results of the perfect foresight model, wherein the optimal land allocation decisions are made based on the values of the crop productivity index in the absence of climate and technology shocks. This deterministic analysis is a useful reference point for further discussion when the uncertainty in food crop yields is introduced. We then present the results of the dynamic stochastic model, where the impact of the intrinsic climate and technology uncertainty is brought into the model optimization stage. Specifically, we generate 1000 sample paths of optimal global land use under different realizations of the stochastic crop productivity index.¹⁰ The results are presented as the difference between the stochastic path and deterministic reference solution.

5.1 Optimal path of global land use under crop yield uncertainty

Figure 3 depicts the optimal allocation of global land use over the next century. The left-hand side of Fig. 3 shows the deterministic paths of different types of land considered in this study, i.e., when the food crop yields are perfectly anticipated. Specifically, it shows three scenarios in which the value of the crop technology index corresponds to (i) expected values of the stochastic crop productivity index (deterministic-baseline scenario), (ii) the most pessimistic climate and bad technology states (deterministic–pessimistic scenario), and (iii) the most optimistic climate and good technology states (deterministic–optimistic scenario). The right-hand side of Fig. 3 shows the difference range between the 1000 simulation paths based on different ex ante realizations of the stochastic crop productivity index and the deterministic-baseline path. The 10 %, 50 %, and 90 % quantile lines represent 10 %, 50 %, and 90 % quantiles of 1000 simulated values, respectively, at each time, and the average line (the + line) represents the average of 1000 simulated values, respectively, at each time.

The right-hand side of Fig. 3 also shows two extreme cases of optimal land use paths conditional on period t realizations of crop productivity index states $A_t(J_{\mathrm{1},t},J_{\mathrm{2},t})$. The realized crop productivity index always takes the highest possible value in a stochastic–optimistic case (the line of upward-pointing triangles) and the lowest possible value in a stochastic–pessimistic case (the line of downward-pointing triangles). As future realizations of the stochastic crop productivity index are uncertain, these extreme stochastic solutions are not the same as the corresponding deterministic solutions, where the values of the future crop productivity index are known with certainty. In the stochastic–optimistic case, for example, potentially lower realizations of future crop productivity index result in larger current-period agricultural land allocation compared to the deterministic–optimistic solution. For other model variables, due to resource limits (e.g., the total land area is unchanged over time) and other constraints, the impact of uncertainty is theoretically difficult to assess. Finally, to ensure consistent interpretation of deterministic and stochastic solutions, the right-hand side of Fig. 3 also shows the difference range between the deterministic–optimistic and deterministic–pessimistic cases relative to their deterministic baselines (using the same markers as on the left-hand side of Fig. 3).

Figure 3 Optimal global land use paths. The solid lines represent the deterministic–baseline solutions. Squares represent the deterministic–optimistic solutions or their distances to the deterministic-baseline solutions. Marks represent the deterministic–pessimistic solutions or their distances to the deterministic-baseline solutions. The shaded areas are the ranges of the sample paths' distances to the deterministic-baseline solutions. Pluses represent the average paths, the dashed lines represent the 10 % quantile, circles represent the median paths, and the dashed–dotted lines represent the 90 % quantile. Upward-pointing triangles and downward-pointing triangles represent the distances of the stochastic–optimistic solutions and the stochastic–pessimistic solutions, respectively, to the deterministic-baseline solutions.

We start a discussion of the left-hand side of Fig. 3, which shows the optimal land use paths under perfect foresight. Given the methodological scope of the paper, we will cover these results briefly. Interested readers should refer to Steinbuks and Hertel (2016) for a more detailed analysis of the perfect foresight model.

Beginning with the description of the baseline scenario, we see that, in the first half of the coming century, the area dedicated to food crops increases, peaking around mid-century and declining significantly thereafter (part 1, top panel) as population growth declines while crop production and food processing technologies improve. Consistent with the literature on 2G biofuel deployment potential (National Research Council, 2011), absent aggressive GHG regulations, and biofuel policies, land allocation for second-generation biofuels remains close to zero until the second half of the coming century (part 1, mid panel). It then becomes viable as fossil fuels become scarce and the costs of producing second-generation biofuels decline. This results in greater land requirements for second-generation biofuel crops. As substitution of pasture land for animal feed in livestock production increases (Taheripour et al., 2013), global pasture area declines while managed forest area increases throughout the entire century (part 1, bottom panel, and part 2, top panel). Finally, unmanaged forest areas decline in response to greater requirements for agricultural land (part 2, mid panel), while protected forest areas will more than double by the end of the coming century in light of strong growth in the demand for ecosystem services (part 2, bottom panel). The other two scenarios exhibit broadly similar dynamics.

We now turn to our main findings about the impacts of uncertainty in the crop productivity index on the distribution of global land resources depicted on the right-hand side of Fig. 3. Compared to deterministic scenarios, this uncertainty results in the additional redistribution of land resources to offset the impact of potentially lower yields. The key reason behind this finding is that social preferences exhibit relative risk aversion (Arrow, 1965; Pratt, 1964) in this stochastic application of the FABLE model. Owing to risk aversion and high adjustment costs of future land conversion, some redistribution takes place as a precautionary policy, even in the absence of actual changes in the states of climate or technology. This is a well-known theoretical result in environmental economics literature (Tsur and Zemel, 2014). Compared to the deterministic baseline scenario, the median (i.e., the 50 % quantile) path of global land use that corresponds to the medium state of climate ($J_{\mathrm{1},t}=\mathrm{3}$) and the medium technological state ($J_{\mathrm{2},t}=\mathrm{2}$) foresees a smaller use of land for food crops (part 1, top panel) and greater use of land for 2G biofuel crops (part 1, mid panel), managed forests (part 2, top panel), and protected land resources (part 2, bottom panel). The differences between stochastic and deterministic paths are small and economically insignificant for all other land resources. This is because land conversion costs of agricultural land for other types of land become larger in the presence of uncertainty. These other types of land have higher adjustment costs of conversion associated with additional time costs of regrowing lumber and livestock and irreversibilities in accessing protected land areas. Land rotation between food crops and 2G biofuel crops is less costly in the FABLE model. This result is consistent with earlier studies that find that closer integration with the energy sector offers greater potential for food–energy substitution and thus also a greater resilience against adverse climate conditions affecting food crop yields (Diffenbaugh et al., 2012; Verma et al., 2014).

While the direction of the effect of the uncertainty in the crop productivity on land conversion can be inferred from the economic theory of environmental and natural resource management under uncertainty (see, e.g., Tsur and Zemel, 2014, and references therein), the extent to which this uncertainty propagates into land conversion depends critically on chosen model structure and parameters. For example, Alexander et al. (2017, p. 1) find that even in the absence of intrinsic uncertainty, “systematic differences in land cover areas are associated with the characteristics modeling approach are at least as great as the differences attributed to scenario variations”. Depending on the assumptions on the substitution of land for other resources, the size of technological progress, and the responsiveness of demand for land-based goods and services to changes in crop productivity, this magnitude can be substantially different from other land use models. However, as we see in Fig. 3, for the same model parameters, the range of land conversion is considerably smaller for the dynamic stochastic model compared to the deterministic scenario analysis. Unlike deterministic analysis, where the solutions paths always stay on a known good or bad trajectory, stochastic analysis allows for future states not to stay in their “best” or “worst” stage. Consequently, optimal land conversion decisions in the stochastic application of the FABLE model are less extreme.

Turning to a specific example, we see from Fig. 3a the difference between the most extreme paths of the stochastic crop productivity index is about 160 million ha by 2100 or about 11 % of the total agricultural area dedicated to food crops. Much of that variation can be attributed to the most extreme (i.e., falling beyond 10th and above 90th percent quantiles) realizations of crop productivity. In line with the argument above, this is because the climate and technological states affecting crop yields are reversible in the stochastic model (that is, if the current state is bad (or good), it could be good (or bad) in the future).

Compared with the deterministic model under the pessimistic (or optimistic) scenario, the social optimum in the stochastic model requires a smaller (or greater) conversion of other types of land to cropland. This is because when the current state of the crop technology index is the worst (best), its future states cannot be worse (better) and have a nonzero probability of being better (worse). The expected future yields will then be better (worse) than the deterministic–pessimistic (optimistic) scenario. As the size of expected crop yields affects the magnitude of the land conversion decisions, the range of stochastic model solutions for agricultural land will be smaller than the range between the most extreme deterministic model solutions. Note this result may not hold if the model solution space is unbounded. This concern does not apply to the stochastic FABLE model because (i) the model’s time horizon is finite, (ii) the crop technology shocks are discrete and finite (hence bounded) based on scientific projections used for the model’s calibration, (iii) all of the model’s state variables (land and fossil fuel resources) are bounded because the total land and the total fossil fuel resources are finite, and (iv) we impose bounds on model decision variables based on the theory of economic dynamics (Barro and Sala-i Martin, 2004), such as strictly positive and finite consumption and output of land-based goods and services; land conversion cannot exceed the total supply of land. The solution space, therefore, must also be bounded because the extent of movement of optimal land uses in any direction is limited by the constraints mentioned above.

Thus, the agricultural land area in the deterministic pessimistic (or optimistic) scenario is larger (or smaller) than the largest (or the smallest) path in the stochastic simulations. For example, in 2100, under the deterministic–pessimistic scenario, the cropland deviation from the deterministic-baseline scenario is about 110 million ha, which is 30 % larger than the largest deviation under the stochastic simulations and more than half of the deviation above the 90 % quantile of the stochastic crop technology index. This result (along with similar results for other land types) demonstrates that scenario analysis can significantly overstate the magnitude of expected agricultural land conversion under uncertain crop yields.

5.2 Optimal path of land-based goods and services under crop yield uncertainty

We now turn to the discussion of optimal paths of land-based goods and services (Fig. 4), keeping the same presentation structure as in the previous section. Deterministic model scenarios (shown on the left-hand side panels of Fig. 4) for land-based goods and services largely mimic trajectories of associated land resources. Specifically, we see steady increases in the production of food crops (part 1, top panel), including livestock and biofuel feedstock, peaking around mid-century and moderating thereafter as consumers satiate their food requirements and the technology of food marketing and processing improves. By 2100, crop production for the livestock feed levels off and even begins to decline. Production of both first- and second-generation biofuels grows as oil becomes more scarce along the baseline path and agricultural yields increase (part 1, mid and bottom panels). Along that optimal path, first-generation biofuels never become a large source of energy consumption, while the production of second-generation biofuels takes off sharply and expands rapidly after 2040 as they become cost-competitive relative to increasingly costly fossil fuels. Production of livestock increases throughout the coming century (part 2, top panel), reflecting shifting diets and the growing demand for processed meat as population income increases (Foley et al., 2011). Production of timber also expands with the growing demand for timber products and further improvements in forest yields (part 2, mid panel). Finally, the consumption of ecosystem services increases throughout most of the coming century as the demand for ecosystem services increases and more natural forest lands become institutionally protected (part 2, bottom panel).

Figure 4 Optimal paths of land-based goods and services. The solid lines represent the deterministic–baseline solutions. Squares represent the deterministic–optimistic solutions or their distances to the deterministic-baseline solutions. Marks represent the deterministic–pessimistic solutions or their distances to the deterministic-baseline solutions. The shaded areas are the ranges of the sample paths' distances to the deterministic-baseline solutions. Pluses represent the average paths, the dashed lines represent the 10 % quantile, circles represent the median paths, and the dashed–dotted lines represent the 90 % quantile. Upward-pointing triangles and downward-pointing triangles represent the distances of the stochastic–optimistic solutions and the stochastic–pessimistic solutions, respectively, to the deterministic-baseline solutions.

The results of the dynamic stochastic model simulations (shown on the right-hand side panels of Fig. 4) show that uncertainty in the crop productivity index has a profound effect on the optimal consumption paths of food crops, biofuels, and livestock and a much smaller effect on consumption paths of merchantable timber and ecosystem services. This result is not very surprising, as crop productivity does not directly affect the production of timber and ecosystem services, whereas indirect land use change effects are relatively small in this stochastic application of the FABLE model. Observe that, unlike land use paths, consumption of land-based goods and services is a decision rather than a state variable, and therefore, there are no significant differences between deterministic and stochastic consumption paths.

Focusing on the larger impacts, we see that between the most extreme paths of the stochastic crop productivity index, the production of food crops varies by about 5.3 billion t along both deterministic and stochastic policy paths. This is a sizable change, which suggests a significant variation in levels of consumption in 2100 along different paths of the stochastic crop productivity index. In the FABLE model, much of the variation in the optimal path of food crops comes on the demand side, with the crop productivity decline resulting primarily in the reduced consumption of processed crops and livestock.¹¹ As shown above, the uncertainty-induced supply response is relatively small along the extensive margin in the dynamic stochastic model (i.e., land conversion). In the Appendix (Fig. A2a), we show that the supply response on the intensive margin is smaller, with the ratio of fertilizers to cropland increasing by fewer than 6 kg ha⁻¹ (or 8 %) under extreme realizations of climate and technology uncertainties. About half of that variation corresponds to the most extreme (i.e., falling beyond 10th and above 90th percent quantiles) realizations of crop productivity. This result indicates that extreme uncertainty in crop productivity could have a significant impact on food consumption over the coming century.

Uncertainty in food crop yields has important implications for the production of the first-generation biofuels that are directly affected by both climate and technology states of food crop yields. The difference between the best and worst states of the crop productivity index is about 31 million tonnes of oil equivalent (Mtoe), which exceeds their expected baseline production in 2100. Although climate and technology states of food crop yields do not directly affect yields of the second-generation biofuel crops, production of second-generation biofuels is nonetheless affected through indirect substitution effects of food for energy in the FABLE demand system. There is a sizable variation in the production of second-generation biofuels between extreme paths of the stochastic crop productivity index, which accounts for 450 Mtoe, or about 30 %, of their total production in 2100.

Given the important contribution of livestock feed in the production of livestock, we can see its production is smaller in the pessimistic scenario and larger in the optimistic scenario. The difference in livestock production between the optimistic and pessimistic scenarios accounts for about 550 million t. As the significance of animal feed in livestock production grows over time, the effect of uncertain crop yields becomes more pronounced. Similar to the result for food crops, the most extreme paths of crop productivity account for about a third of all variation in livestock production.

6 Conclusions

This paper shows the effects of uncertainties associated with nonstationary biophysical processes and technological change on the optimal allocation of natural resources in the long run. In doing so, it applies SCEQ, a cutting-edge computational method for solving nonstationary dynamic high-dimensional stochastic problems, to FABLE, a multi-sectoral dynamic model of global land use.

The study focuses on uncertainty in future crop yields, one of the core uncertainties affecting the evolution of global land use in the long run. Combining scenarios from global climate models and high-resolution output from spatial crop-simulation models for four major crops, it comes up with a plausible range of realizations of climate shocks and their effect on future crop yields. These estimates are supplemented with an extensive survey of recent agro-economic and biophysical studies assessing the potential for closing yield gaps, as well as attaining further advances in potential yields through plant breeding.

The paper's key insight is to illustrate the magnitude of optimal land conversion decisions in the context of different realizations of stochastic crop productivity. Consistent with the economic theory of natural resource management under uncertainty, the agricultural productivity shocks from either adverse climate impacts or unexpected limits on further technological progress result in additional conversion of scarce land resources to offset the impact of potentially lower yields. Owing to intertemporal substitution, some of that conversion takes place even in the absence of actual realization of the climate shocks or technology outcomes. This expansion is accompanied by changes in the consumption of processed food, livestock, and biofuel- and land-based products most affected by changes in crop productivity.

The chosen model (FABLE) seeks to balance computational complexity and economic tractability. Similar to other dynamic economic models, it assumes a bounded solution space, excluding unrealistic scenarios of infinitely low or high crop productivity. It also ignores many features that are standard in more advanced computational land and other resource-use models. Future research should focus on integrating economic decisions under uncertainty into large dynamic natural resource models that feature spatial disaggregation at the regional or zonal level, a more extensive representation of the energy sector, and different types of resources and their production derivatives. Another promising research direction would be to incorporate a more detailed representation of uncertain states backed by an econometric analysis that recovers underlying distributions of uncertain natural resource drivers over time.

Appendix A: FABLE model description

This section describes key elements of the FABLE model, as well as its equations, variables, and model parameters. For a full description of the model, including details on model baseline calibration and extensive sensitivity analysis, please refer to Hertel et al. (2016) and Steinbuks and Hertel (2016) and technical appendices therein.

A1 Primary resources

Primary resources comprise land, liquid fossil fuels, and other primary inputs, e.g., labor and capital. The supply of land is fixed and faces competing uses that are determined endogenously by the model. The flow of liquid fossil fuels evolves endogenously along their optimal path, accounting for exogenous discoveries in new fossil fuel reserves. The endowment of other primary inputs is exogenous and evolves along the prespecified global economy growth path.

A1.1 Land

The total land endowment in the model, L^total, is fixed. It belongs to a global planner, which fully redistributes land rents back to consumers of land use goods and services. For each period of time t there are four profiles of land in the economy. They include unmanaged forest land, L^N, agricultural land, L^A, pasture land, L^P, and commercially managed forest land, L^C. The agricultural land area can be allocated for the cultivation of food crops (denoted L^A,c), and second-generation biofuel feedstocks (denoted L^A,b2). We assume that the natural forest land consists of two types. Institutionally protected land, L^R, includes natural parks, biodiversity reserves, and other types of protected forests. This land is used to produce ecosystem services for society and cannot be converted to commercial land. Unmanaged natural land, L^N, can be accessed and either converted to managed land or to protected natural land. Once the natural land is converted to managed land, its potential to yield ecosystem services is diminished. This potential can be partially restored for managed forests with significant land rehabilitation costs incurred. The use of managed land can be shifted between cropland, forestland, and pasture land (see Fig. 1 for a graphical representation of these transitions). We denote land transition flows from land type i to land type j as Δ^i,j (a negative value means a transition from land type j to land type i). Equations describing the allocation of land across time and different uses are as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{L}^{\mathrm{total}}=\sum _{i=\mathrm{A},\mathrm{P},\mathrm{C},\mathrm{N},\mathrm{R}}{L}_t^i,\text{(A2)}& {\displaystyle }{L}^{\mathrm{A}}=L^{\mathrm{A},\mathrm{c}}+L^{\mathrm{A},\mathrm{b}\mathrm{2}},\text{(A3)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{N}}=L_t^{\mathrm{N}}-{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}-{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}+{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}},\text{(A4)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{A}}=L_t^{\mathrm{A}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}-{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}}+{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}},\text{(A5)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{P}}=L_t^{\mathrm{P}}+{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}},\text{(A6)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{R}}=L_t^{\mathrm{R}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}.\end{array}$$

Equations (A1) and (A2) define, respectively, the composition of total land and agricultural land in the economy. Equations (A3)–(A5) describe the transitions for unmanaged land, agricultural land, and pasture land.¹² Equation (A6) shows the growth path of protected natural land.

Accessing the natural lands comes at a cost associated with building roads and other infrastructure (Golub et al., 2009). In addition, converting natural land to reserved land entails additional costs associated with passing legislation to create new natural parks. We denote the natural land access, rehabilitation, and protection costs as $C^{\mathrm{N},\mathrm{A},\mathrm{R}},$ $C^{\mathrm{C},\mathrm{N}},$ and $C^{\mathrm{N},\mathrm{R}},$, respectively. There are also costs of switching between the cropland and the pasture land, denoted as C^A,P. We assume that all these costs are continuous, monotonically increasing, and strictly convex functions of converted land. Since we are not aware of empirical studies estimating the magnitudes of long-term adjustment costs in land conversion problems, we choose to calibrate these parameters to match historical land conversion patterns. There are no additional costs of natural land conversion to commercial land, as the revenues from deforestation offset these costs.

Managed forests are characterized by v_max vintages of tree species with vintage ages $v=\mathrm{1},\mathrm{\dots },v_{max}$. At the end of period t, each hectare of managed forest land, $L_{v,t}^{\mathrm{C}},$ has an average density of tree vintage age v, with the initial allocation given and denoted by $L_{v,\mathrm{0}}^{\mathrm{C}}$. The forest rotation ages and management are endogenously determined. Each period the managed forest land can be either planted, harvested, or left to mature. The newly planted trees occupy Δ^C,C hectares of land and reach the average age of the first tree vintage next period. The harvested area of tree vintage age v occupies $\mathrm{\Delta }{}_v^{\mathrm{C},\mathrm{H}}$ hectares of forest land. The difference between the harvested area of all tree vintage ages and the newly planted area is used for cropland, i.e.,

$${\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}}=\sum _v{\mathrm{\Delta }}_{v,t}^{\mathrm{C},\mathrm{H}}-{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}}.$$

The following equations describe land use of managed forests:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{L}_t^{\mathrm{C}}=\sum _{v=\mathrm{1}}^{v_{max}}{L}_{v,t}^{\mathrm{C}},\text{(A8)}& {\displaystyle }{L}_{v+\mathrm{1},t+\mathrm{1}}^{\mathrm{C}}=L_{v,t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v,t}^{\mathrm{C},\mathrm{H}},v<v_{max}-\mathrm{1},\text{(A9)}& {\displaystyle }\begin{array}{rl}{L}_{v_{max},t+\mathrm{1}}^{\mathrm{C}}& =L_{v_{max},t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v_{max},t}^{\mathrm{C},\mathrm{H}}-{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}}\\ & +L_{v_{max}-\mathrm{1},t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v_{max}-\mathrm{1},t}^{\mathrm{C},\mathrm{H}},\end{array}\text{(A10)}& {\displaystyle }{L}_{\mathrm{1},t+\mathrm{1}}^{\mathrm{C}}={\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}}.\end{array}$$

Equation (A7) describes the composition of managed forest area across vintages. Equation (A8) illustrates the harvesting dynamics of forest areas with the ages $v_{max}-\mathrm{1}$ and v_max. Equation (A10) shows the transition from the planted area to new forest vintage area.

The average harvesting and planting costs per hectare of new forest planted, c^o,H and $c^{\mathrm{o},\mathrm{C}},$ are invariant to scale and are the same across all vintages. Harvesting managed forests and conversion of harvested forest land to agricultural land is subject to additional near-term adjustment costs, c^H. The specific functional forms of land conversion costs are shown in Appendix C and Eqs. (C32)–(C37).

Thus, we have defined the vector of land state variables,

$$\mathbf{L}=\left(L^{\mathrm{N}},L^{\mathrm{A}},L^{\mathrm{P}},L^{\mathrm{R}},L_{\mathrm{1}}^{\mathrm{C}},\mathrm{\dots },L_{v_{max}}^{\mathrm{C}}\right),$$

and its associated transition laws.

A1.2 Fossil fuels

The initial stock of liquid fossil fuels, X^F, is exogenous, and each period of time t adds a new number of fossil fuels, ${\mathrm{\Delta }}^{\mathrm{F},\mathrm{D}},$ which reflects exogenous technological progress in fossil fuel exploration. This technological progress comprises of both discoveries on new exploitable oil and gas fields, as well as development of new technologies for extraction of non-conventional fossil fuels. The economy extracts fossil fuels, which have two competing uses in our partial-equilibrium model of land use. A part of the extracted fossil fuels, ${\mathrm{\Delta }}_t^{\mathrm{F},n}$, is converted to fertilizers that are further used in the agricultural sector. The remaining number of fossil fuels, ${\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}},$ is combusted to satisfy the demand for energy services. The following equation describes supply of fossil fuels:

$$\begin{array}{}\text{(A11)}& X_{t+\mathrm{1}}^{\mathrm{F}}=X_t^{\mathrm{F}}-{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}-{\mathrm{\Delta }}_t^{\mathrm{F},n}+{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{D}}.\end{array}$$

The cost of fossil fuels, c^F, reflects the expenditures for fossil fuel extraction, refining, transportation, and distribution, as well the costs associated with emissions control (e.g., Pigouvian taxes) in the non-land-based economy. We assume that the cost of fossil fuels is a nonlinear quadratic function with accelerating costs as the stock of fossil fuels depletes (Nordhaus and Boyer, 2000):

$$\begin{array}{}\text{(A12)}& c_t^{\mathrm{F}}={\mathit{\xi }}_{\mathrm{1}}^{\mathrm{F}}{\left({\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}+{\mathrm{\Delta }}_t^{\mathrm{F},n}\right)}^{\mathrm{2}}\left({\displaystyle \frac{X_{\mathrm{0}}^{\mathrm{F}}+{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{D}}}{X_t^{\mathrm{F}}+{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{D}}}}\right),\end{array}$$

where the parameter ${\mathit{\xi }}_{\mathrm{1}}^{\mathrm{F}}$ captures the curvature of the liquid fossil fuel cost function.

A1.3 Other primary resources

The initial endowment of all other primary resources in the non-land-based economy, such as labor, physical and human capital, and material inputs, X^O, is exogenous in this model. We assume that the growth rate of all other primary resources is a weighted average of the population growth, which reflects demographic changes, and the physical capital growth, κ^o,X. The following equation describes the supply of other primary inputs:

$$\begin{array}{}\text{(A13)}& X_t^{\mathrm{O}}=X_{\mathrm{0}}^{\mathrm{O}}\left[{\mathit{\alpha }}^{\mathrm{o},\mathrm{l}}{\displaystyle \frac{{\mathrm{\Pi }}_t}{{\mathrm{\Pi }}_{\mathrm{0}}}}+\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{o},\mathrm{l}}\right){\left(\mathrm{1}+{\mathit{\kappa }}^{o,X}\right)}^t\right],\end{array}$$

where Π_t is the economy's population, and α^o,l is the share of population growth to the growth rate of all other primary resources. Other primary inputs can be used for the production of land-based goods and services or be converted to final goods and services in the non-land economy. The production costs incurred using these inputs are exogenous and have an “iceberg” representation; i.e., they are subtracted from the gross output of land-based goods and services. Thus, state variables for resources other than land are defined as

$$\mathbf{X}=(X^{\mathrm{F}},X^{\mathrm{O}}).$$

As X^O is exogenous and deterministic, it is a degenerated state variable and not counted as a state variable for model solution purposes.

A2 Intermediate inputs

We analyze six intermediate inputs used in the production of land-based goods and services: petroleum products, fertilizers, crops, biofuels, and raw timber. Fossil fuels are refined and converted to either petroleum products, x^p, that are further combusted or to fertilizers, xⁿ, that are used to boost yields in the agricultural sector. Agricultural land and fertilizers are combined to grow food crops, x^c or 2G biofuel crops, x^c,b2. Food crops can be further converted into processed food and 1G biofuels, x^b1, or used as an animal feed, x^c,l. 2G biofuel crops can only be converted into 2G biofuels, x^b2. 1G biofuels substitute imperfectly for liquid fossil fuels in final energy demand, whereas 2G biofuels and liquid fossil fuels are the perfect substitutes The food crops used as animal feed and pasture land are combined to produce raw livestock, x^l. Harvesting managed forests yields raw timber, x^w, that is further used in timber processing. The production functions for intermediate inputs can be illustrated by the following equations

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{x}_t^j& =g^j\left({\mathrm{\Delta }}_t^{\mathrm{F},\mathit{\{}\mathrm{E},n\mathit{\}}},L_t^{\{\mathrm{A},\mathrm{P}\}},\sum _v\mathrm{\Delta }{}_{v,t}^{\mathrm{C},\mathrm{H}},x^{\mathrm{c},\{\mathrm{l},\mathrm{b}\}}\right),\\ & j=p,n,c,b,l,w,\end{array}\end{array}$$

where ${\mathrm{\Delta }}_t^{\mathrm{F},\mathit{\{}\mathrm{E},n\mathit{\}}}$ represents that either ${\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}$ or ${\mathrm{\Delta }}_t^{\mathrm{F},n}$ is an argument of g^j, similarly for $L_t^{\{\mathrm{A},\mathrm{P}\}}$ and $x^{\mathrm{c},\{\mathrm{l},\mathrm{b}\}}$. The specific functional forms of g^j(⋅) are shown in Appendix C and Eqs. (C13)–(C21).

A3 Final goods and services

We consider five per capita land-based services that are consumed in the final demand: services from processed crops, y^f, livestock, y^l, energy, y^e, timber, y^w, and ecosystem services, y^r. Processed crops, livestock, and timber are, respectively, products of food crops, raw livestock, and timber processing. The production of energy services combines liquid fossil fuels with biofuels, and the resulting mix is further combusted. The ecosystem services are the public good to society, which captures recreation, biodiversity, and other environmental goods and services. To close the demand system, we also include other goods and services, y^o, which comprise the consumption of other primary inputs not spent on the production of land-based goods and services. We have defined all state variables for the deterministic model:

$$\mathit{S}:=(\mathbf{L},\mathbf{X}),$$

and the vector of decision variables

$$\begin{array}{rl}{\mathit{a}}_t:& =({\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}},{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}},{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}},{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}},{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}},{\mathrm{\Delta }}_{\mathrm{1},t}^{\mathrm{C},\mathrm{H}},\mathrm{\dots },\\ & {\mathrm{\Delta }}_{v_{max},t}^{\mathrm{C},\mathrm{H}},{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}},{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}},{\mathrm{\Delta }}_t^{\mathrm{F},n},L_t^{\mathrm{A},\mathrm{F}},L_t^{\mathrm{A},\mathrm{B}},{\mathit{x}}_t,{\mathit{y}}_t),\end{array}$$

where ${\mathit{x}}_t\equiv [x_t^{\mathrm{p}},x_t^n,x_t^{\mathrm{c}},x_t^{\mathrm{b}},x_t^{\mathrm{l}},x_t^{\mathrm{w}},x_t^{\mathrm{c},\mathrm{l}},x_t^{\mathrm{c},\mathrm{b}}]$ and ${\mathit{y}}_t=\left(y_t^{\mathrm{f}},y_t^{\mathrm{l}},y_t^{\mathrm{e}},y_t^{\mathrm{w}},y_t^{\mathrm{r}},y_t^{\mathrm{o}}\right)$. The production functions for final per capita land-based goods and services can be illustrated by the following equation:

$$\begin{array}{}\text{(A15)}& y_t^i={\mathcal{Y}}_t^i\left({\mathit{S}}_t,{\mathit{a}}_t\right),i=f,l,e,w,r,\end{array}$$

where some arguments in ${\mathcal{Y}}_t^i\left(\cdot \right)$ could be redundant. It follows from Eq. (A15) that the production of final goods and services involves the combination of land resources and intermediate inputs. The specific functional forms of ${\mathcal{Y}}_t^i\left(\cdot \right)$ are shown in Appendix C and Eqs. (C22)–(C27), which are functions of L and {x^j}. All these equations constitute a part of the feasibility constraint a_t∈𝒟_t(S_t).

The production of intermediate inputs or final land-based goods and services i incurs costs, c^o,i, that are subtracted from other available primary resources. The remaining number of other primary resources is converted into other goods and services, which are subsequently consumed in the final demand. As the focus of this model is on the utilization of land-based resources, we introduce the other goods and services, y^o, in a very simplified manner. We introduce no additional cost of producing other goods and services, assuming that it is reflected in the size of the endowment of other primary inputs. The specific functional form for y^o is shown in Appendix C, Eq. (C26).

A4 Preferences

The economy's per capita utility, u, is derived from the per capita consumption of processed crops, livestock, timber, energy and ecosystem services, and other goods and services. Following the macroeconomic literature, we assume constant relative risk-aversion utility,

$$\begin{array}{}\text{(A16)}& u\left(\mathit{y}\right)={\displaystyle \frac{\mathcal{C}(\mathit{y}{)}^{\mathrm{1}-\mathit{\gamma }}}{\mathrm{1}-\mathit{\gamma }}},\end{array}$$

where y is the per capita consumption bundle of goods and services, 𝒞(y) is a nonlinear aggregator over y, and γ is the coefficient of relative risk aversion, which captures the economy's attitude to uncertain events. We choose a non-homothetic AIDADS preference (Rimmer and Powell, 1996) to compute 𝒞(y) implicitly:

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}\log \left(\mathcal{C}\left(\mathit{y}\right)\right)& =\sum _{q=f,l,e,w,r,o}\left({\displaystyle \frac{{\mathit{\alpha }}_q+{\mathit{\beta }}_q\mathcal{C}\left(\mathit{y}\right)}{\mathrm{1}+\mathcal{C}\left(\mathit{y}\right)}}\right)\\ & \log \left(y^q-{\underset{\mathrm{¯}}{y}}^q\right),\end{array}\end{array}$$

where α, β, and ${\underset{\mathrm{¯}}{y}}^q$ are positive parameters with ${\sum }_q{\mathit{\alpha }}_q={\sum }_q{\mathit{\beta }}_q=\mathrm{1}$. These preferences place greater value on ecosystem services and smaller value on additional consumption of food, energy, and timber products as society becomes wealthier. When γ=1, our utility function is equivalent to the AIDADS utility.

A5 Welfare

We denote the transition laws of land, in Eqs. (A3)–(A6) and Eqs. (A8)–(A10), as

$$\begin{array}{}\text{(A18)}& {\mathbf{L}}_{t+\mathrm{1}}={\mathit{G}}_t^{\mathrm{L}}({\mathit{S}}_t,{\mathit{a}}_t),\end{array}$$

and the transition laws for other resources, in Eqs. (A11)–(A13), as

$$\begin{array}{}\text{(A19)}& {\mathbf{X}}_{t+\mathrm{1}}={\mathit{G}}_t^X({\mathit{S}}_t,{\mathit{a}}_t).\end{array}$$

Combining Eqs. (A18) and (A19), we have

$$\begin{array}{}\text{(A20)}& {\mathit{S}}_{t+\mathrm{1}}={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t)\end{array}$$

for the deterministic model in the notations of Sect. 2.

The objective of the planner is to maximize the total expected welfare, which is the cumulative expected utility (i.e., a sum of utilities in each state weighted by the probability of each state in a given period) of the population's consumption of final goods and services, y, discounted at the constant rate $\mathit{\delta }>\mathrm{0}.$ The planner allocates managed agricultural, pasture, and forest lands for crop, livestock, and timber production, the scarce fossil fuels, and protected natural forests to solve the following problem:

$$\begin{array}{}\text{(A21)}& \underset{\mathit{a}}{max}\sum _{t=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\delta }}^t\mathcal{U}({\mathit{S}}_t,{\mathit{a}}_t),\end{array}$$

subject to the transition laws (A20) and the feasibility constraints,

$${\mathbf{F}}_t({\mathit{S}}_t,{\mathit{a}}_t),\ge \mathrm{0}$$

which include Eqs. (A15), (C22)–(C27), (C26), and (A17) and nonnegativity constraints for the variables. Here,

$$\begin{array}{}\text{(A22)}& \mathcal{U}({\mathit{S}}_t,{\mathit{a}}_t)=u\left({\mathit{y}}_t\right){\mathrm{\Pi }}_t\end{array}$$

is the utility function in the notations of Sect. 4.

Appendix B: Quantifying the uncertainty in crop yields

B1 Uncertainty in agricultural technology

Advances in crop technology are very difficult to predict due to four interconnected factors (Fischer et al., 2011). First, there is significant uncertainty about the potential for exploiting large and economically significant yield gaps (i.e., the differences between observed and potential crop yields) in developing countries, especially those in sub-Saharan Africa. A second and closely related point is that it is unclear how fast available yield-enhancing technologies can be adopted at a global scale. Third, there is a significant variation in developing countries' institutions and policies that make markets work better and provide a conducive environment for agricultural technology adoption. Finally, while plant breeders continue to make steady gains in further advancing crop yields, progress depends on the level of funding provided for agricultural research. This has proven to be somewhat volatile, with per capita funding falling in the decades leading up to the recent food crisis (Alston and Pardey, 2014). Food prices have risen since 2007, which has stimulated new investments. However, whether this interest will be sustained remains to be seen. Overall, progress from conventional breeding is becoming more difficult. Transgenic (genetic modification) technologies have a proven record of more than a decade of safe and environmentally sound use and thus offer huge potential to address critical biotic and abiotic stresses in the developing world. However, expected yield gains, costs of further developing these technologies, and the political acceptance of genetically modified foods are all highly uncertain.

To quantify the extent to which the advances in crop technology can further boost agricultural yields over the next century, we first need to assess the magnitude of existing yield gaps at the global scale. In a comprehensive study, Lobell et al. (2009) report a significant variation in the ratios of actual to potential yields for major food crops across the world, ranging from 0.16 for tropical lowland maize in sub-Saharan Africa to 0.95 for wheat in Haryana, India. For the purposes of this study, we employ the results of Licker et al. (2010), who conduct comprehensive yield gap analysis using a global crop dataset of harvested areas and yields for 175 crops on a 0.5° geographic grid of the planet for the year 2000. Using these estimates, we calculate the global yield gap as the grid-level output-weighted yield gap of the four most important food crops (wheat, maize, soybeans, and rice). The resulting estimate suggests that average yields are 53 % of potential yields, which is close to the median estimates by Lobell et al. (2009). As a further robustness check, we employ the Decision Support System for Agrotechnology Transfer (DSSAT) crop-simulation model (Jones et al., 2003), run globally on a 0.5° grid in the parallel System for Integrating Impacts Models and Sectors (pSIMS; Elliott et al., 2014b) to simulate yields of the same four major food crops under the best agricultural management conditions and compare simulated yields to their observed yields. The resulting yield gap estimates were not substantially different.

In the optimistic (i.e., good) state of advances in crop technology, we assume that yields will continue to grow linearly throughout the coming century, eliminating the yield gap by 2100. This high-yield scenario rests on the assumption of continued strong growth in investment in agricultural research and development, widespread acceptance of genetically modified crops, continuing institutional reforms in developing countries, and public and private investments in the dissemination of new technologies. The erosion of any one of these component assumptions will likely result in a slowing of crop technology improvements. And there are some grounds for pessimism. In a comprehensive statistical analysis of historical crop production trends, Grassini et al. (2013) note that

despite the increase in investment in agricultural R&D and education … the relative rate of yield gain for the major food crops has decreased over time together with evidence of upper yield plateaus in some of the most productive domains. For example, investment in R&D in agriculture in China has increased threefold from 1981 to 2000. However, rates of increase in crop yields in China have remained constant in wheat, decreased by 64 % in maize as a relative rate and are negligible in rice. Likewise, despite a 58 % increase in investment in agricultural R&D in the United States from 1981 to 2000 (sum of public and private sectors), the rate of maize yield gain has remained strongly linear.

To capture the possibility of much slower technological improvement in the coming century, we specify two more pessimistic scenarios. In the medium state of technology, rather than closing the yield gap by 2100, average yields in 2100 are just three-quarters of yield potential at that point in time. In the bad state of technology, there is no technological progress, and the crop yields stay the same as at the beginning of the coming century. This is the path on which we begin the simulation in 2004. As previously noted, we then specify probabilities with which the crop technology index evolves across the different states of technology.

B2 Uncertainty in climate change impacts

In addition to crop technology uncertainty, there is great uncertainty about the physical environment in which this technology will be deployed. In particular, long-run changes in temperature and precipitation are likely to have an important impact on land productivity in agriculture (IPCC, 2014) and, therefore, the global pattern of land use. Quantification of the impact of climate change on agricultural yields requires coming to grips with three interconnected factors (Alexandratos, 2011). First, there is significant uncertainty in future GHG concentrations along the long-run growth path of the global economy. Second, the general circulation models (GCMs) developed by climate scientists to translate these uncertain GHG concentrations into climate outcomes disagree about the spatially disaggregated deviations of temperature and precipitation from baseline levels. Finally, there is significant uncertainty in the biophysical models used to determine how changes in temperature and precipitation will affect plant growth and the productivity of agriculture in different agroecological conditions. The impact of climate change on food crop yields depends critically on their phenological development, which, in turn, depends on the accumulation of heat units, typically measured as growing degree days (GDDs). More rapid accumulation of GDDs due to climate change speeds up phenological development, thereby shortening key growth stages, such as the grain-filling stage, hence reducing potential yields (Long, 1991). However, rising concentrations of CO₂ in the atmosphere result in an increase in potential yields due to improved water use efficiency, often dubbed the “CO₂ fertilization effect” (Long et al., 2006). Sorting out the relative importance of these effects and achieving greater confidence in evaluations of climate impacts on agricultural yields remains an important research question in the agronomic literature (Cassman et al., 2010; Rosenzweig et al., 2014).

To quantify the uncertainty in climate impacts on agricultural yields, we use ISIMIP fast-track crop-simulation model data (Rosenzweig et al., 2014). Specifically, we obtain results of four crop-simulation models: GEPIC (Liu et al., 2007), LPJmL (Bondeau et al., 2007), pDSSAT (Jones et al., 2003), and PEGASUS (Deryng et al., 2011). All models are run globally on a 0.5° grid over the period between 1971 and 2099 and weighted by the output of four major food crops (maize, soybeans, wheat, and rice). Since the FABLE model is not spatially explicit, we further aggregate gridded crop yields at a global scale using weights of the size of the aggregate crop output per grid cell. To ensure the simulation results' comparability with the structural parameters of the FABLE model, all models are run under the Representative Concentration Pathway 6.0 W m⁻² (RCP6) GHG forcing scenario (Moss et al., 2008). We also consider alternative assumptions on CO₂ fertilization effects. Observe that our results are based on four crop-simulation models, though Rosenzweig et al. (2014) consider seven crop-simulation models. The remaining three models have fewer crops and/or temporal frames for model baseline and are thus omitted. Rosenzweig et al. (2014) find that five models, including GEPIC, LPJmL, and pDSSAT models considered in this analysis, yield broadly similar predictions. One model (LPJ–GUESS) not covered here has much higher variation in predicted crop yields under different climate scenarios. Our results may, therefore, understate the range of uncertainty in the climate change impacts on potential crop yields.

To quantify uncertainty in temperature increases due to climate change, we employ outputs for five global climate models (GCMs): GFDL-ESM2M (Dunne et al., 2013), HadGEM2-ES (Collins et al., 2008), IPSL-CM5A-LR (Dufresne et al., 2013), MIROC-ESM-CHEM (Watanabe et al., 2011), and NorESM1-M (Bentsen et al., 2013). For each of the simulations, we fit a linear trend to parsimoniously characterize the evolution of crop yields in the face of climate change over the coming century.

Figure C1 summarizes simulation results for four crop-simulation models and five climate models (with and without fertilization effects) in 2100, normalized relative to assumed yield potential in the absence of climate change. There is significant heterogeneity in terms of both direction and magnitude of climate impacts on agricultural yields across global climate models when the CO₂ fertilization effect is considered.¹³ Regardless of the chosen climate model, for the scenario with fertilization effects, two out of four crop-simulation models (LPJmL and pDSSAT) predict a moderate increase in potential yields (5 %–15 %), whereas the PEGASUS model predicts a large decline in potential yields (20 %–30 %). The GEPIC model predicts that on average crop yields will be little changed, showing a small increase in crop yields for some climate models and a small decline for other models. The predictions of LPJmL and pDSSAT models are reversed when CO₂ fertilization effects are removed, showing a decline of about 10 %–15 % in potential yields. The PEGASUS model predicts an even larger decline in potential yields (30 %–35 %), whereas the predictions of GEPIC model show a moderate decline of about 5 %–10 % in potential yields.

Given a large variation in model predictions, we construct five states for potential crop yields under uncertain climate change. These states correspond to quintiles of the distribution of different model outcomes for potential crop yields by 2100. Under two optimistic states of the world, we observe a 2 % and 15 % increases in potential crop yields relative to the model baseline, whereby significant CO₂ fertilization effects offset the negative effects of climate change. For the next two states, we see 15 % and 19 % declines in potential crop yields relative to model baseline, whereby CO₂ fertilization effects are either small or nonexistent, and the negative effects of climate change tend to prevail. Finally, under the most pessimistic states of the world, drastic adverse effects of climate change combined with the absence of any CO₂ fertilization effects result in a 36 % decline in potential crop yields relative to the model baseline.

B3 Transition probabilities

The five possible values of the climate state J_1,t are ${\mathcal{J}}_{\mathrm{1},\mathrm{1}}=\mathrm{0.64}$, ${\mathcal{J}}_{\mathrm{1},\mathrm{2}}=\mathrm{0.85}$, ${\mathcal{J}}_{\mathrm{1},\mathrm{3}}=\mathrm{0.89}$, ${\mathcal{J}}_{\mathrm{1},\mathrm{4}}=\mathrm{1.02}$, and ${\mathcal{J}}_{\mathrm{1},\mathrm{5}}=\mathrm{1.15}$, and its probability transition matrix is

$$P_{\mathrm{1}}=\left[\begin{array}{ccccc}\mathrm{0.5}& \mathrm{0.25}& & & \\ \mathrm{0.5}& \mathrm{0.5}& \mathrm{0.25}& & \\ & \mathrm{0.25}& \mathrm{0.5}& \mathrm{0.25}& \\ & & \mathrm{0.25}& \mathrm{0.5}& \mathrm{0.5}\\ & & & \mathrm{0.25}& \mathrm{0.5}\end{array}\right],$$

where $P_{\mathrm{1},i,j}$ represents the probability from the jth value of J_1,t to the ith value for $\mathrm{1}\le i,j\le \mathrm{5}$. The three possible values of the technological state J_2,t are ${\mathcal{J}}_{\mathrm{2},\mathrm{1}}=\mathrm{1.45}$, ${\mathcal{J}}_{\mathrm{2},\mathrm{2}}=\mathrm{1.675}$, and ${\mathcal{J}}_{\mathrm{2},\mathrm{3}}=\mathrm{1.9}$, and its probability transition matrix is

$$P_{\mathrm{2}}=\left[\begin{array}{ccc}\mathrm{0.4423}& \mathrm{0.1416}& \mathrm{0.1311}\\ \mathrm{0.4139}& \mathrm{0.669}& \mathrm{0.4367}\\ \mathrm{0.1438}& \mathrm{0.1894}& \mathrm{0.4322}\end{array}\right],$$

where $P_{\mathrm{2},i,j}$ represents the probability from the jth value of J_2,t to the ith value for $\mathrm{1}\le i,j\le \mathrm{3}$. We assume that J_2,t is independent of J_1,t.

B4 Model

After we add the risks, the state vector becomes

$$\mathit{S}:=(\mathbf{L},\mathbf{X},\mathbf{J}),$$

where ${\mathbf{J}}_t=(J_{\mathrm{1},t},J_{\mathrm{2},t})$. And J is a Markov chain, so it can be represented as ${\mathbf{J}}_{t+\mathrm{1}}={\mathit{G}}_t^J({\mathbf{J}}_t,{\mathit{ϵ}}_t)$, where ϵ_t is a vector of shocks with zero means. The problem is

$$\begin{array}{}\text{(B1)}& \underset{\mathit{a}}{max}\mathbb{E}\left\{\sum _{t=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\delta }}^t\mathcal{U}({\mathit{S}}_t,{\mathit{a}}_t)\right\},\end{array}$$

subject to

$$\begin{array}{ll}{\displaystyle }{\mathbf{L}}_{t+\mathrm{1}}& {\displaystyle }={\mathit{G}}_t^{\mathrm{L}}({\mathit{S}}_t,{\mathit{a}}_t),\\ {\displaystyle }{\mathbf{X}}_{t+\mathrm{1}}& {\displaystyle }={\mathit{G}}_t^X({\mathit{S}}_t,{\mathit{a}}_t),\\ {\displaystyle }{\mathbf{J}}_{t+\mathrm{1}}& {\displaystyle }={\mathit{G}}_t^J({\mathbf{J}}_t,{\mathit{ϵ}}_t),\end{array}$$

and a_t∈𝒟_t(S_t) representing the feasibility constraints, that is, inequality constraints and the equations other than the above transition laws. The above transition laws are just a special case of

$${\displaystyle }{\mathit{S}}_{t+\mathrm{1}}{\displaystyle }={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t,{\mathit{ϵ}}_t)$$

in the notations of Sect. A2 of this Appendix, so we can implement the SCEQ method to solve the dynamic stochastic programming problem. Since our time of interest is $T^{*}=\mathrm{100}$ years, we change the problem (B1) to have a finite horizon with T=400 years, as a larger T has little impact on our solution in the first 100 years.

In the step 2 of Algorithm 1 for the solution at time s, we replace ϵ_t by its zero mean to have ${\mathit{S}}_{t+\mathrm{1}}={\mathit{G}}_t({\mathit{S}}_t,{\mathit{a}}_t,\mathrm{0})$; that is, ${\mathbf{J}}_{t+\mathrm{1}}={\mathit{G}}_t^J({\mathbf{J}}_t,\mathrm{0})$. But this ${\mathbf{J}}_{t+\mathrm{1}}={\mathit{G}}_t^J({\mathbf{J}}_t,\mathrm{0})$ is only for simplicity in notations. In fact, since J is a Markov chain, we replace J_t by its mean conditional on the realized value of J_s (i.e., its certainty equivalent approximation):

$$\left[{\mathcal{J}}_{\mathrm{1}}{\mathit{\pi }}_{\mathrm{1},t,s},{\mathcal{J}}_{\mathrm{2}}{\mathit{\pi }}_{\mathrm{2},t,s}\right]$$

for all t≥s, where ${\mathcal{J}}_{\mathrm{1}}=({\mathcal{J}}_{\mathrm{1},\mathrm{1}},\mathrm{\dots },{\mathcal{J}}_{\mathrm{1},\mathrm{5}})$, ${\mathcal{J}}_{\mathrm{2}}=({\mathcal{J}}_{\mathrm{2},\mathrm{1}},{\mathcal{J}}_{\mathrm{2},\mathrm{2}},{\mathcal{J}}_{\mathrm{2},\mathrm{3}})$, ${\mathit{\pi }}_{\mathrm{1},t,s}$, and ${\mathit{\pi }}_{\mathrm{2},t,s}$ are two-column vectors representing probability distributions of J_1,t and J_2,t conditional on the realized values of J_1,s and J_2,s, respectively. If the realized values of J_1,s and J_2,s are 𝒥_1,i and 𝒥_2,j, respectively, then we have ${\mathit{\pi }}_{\mathrm{1},t,s}=P_{\mathrm{1}}^{t-s}{\mathit{\pi }}_{\mathrm{1},s,s}$ and ${\mathit{\pi }}_{\mathrm{2},t,s}=P_{\mathrm{2}}^{t-s}{\mathit{\pi }}_{\mathrm{2},s,s}$, where ${\mathit{\pi }}_{\mathrm{1},s,s}$ is a length 5 column vector, with 1 at the ith element and 0 everywhere else, and ${\mathit{\pi }}_{\mathrm{2},s,s}$ is a length 3 column vector, with 1 at the jth element and 0 everywhere else.

Appendix C: Model equations, variables and parameters

C1 Equations

C1.1 Land use

$$\begin{array}{}\text{(C1)}& {\displaystyle }L=\sum _{i=\mathrm{A},\mathrm{P},\mathrm{C},\mathrm{N},\mathrm{R}}{L}_t^i\text{(C2)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{N}}=L_t^{\mathrm{N}}-{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}-{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}+{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}}\text{(C3)}& {\displaystyle }{L}_t^{\mathrm{A}}=L_t^{\mathrm{A},\mathrm{c}}+L_t^{\mathrm{A},\mathrm{b}\mathrm{2}}\text{(C4)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{A}}=L_t^{\mathrm{A}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}-{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}}+{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}}\text{(C5)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{P}}=L_t^{\mathrm{P}}+{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}}\text{(C6)}& {\displaystyle }{L}_{t+\mathrm{1}}^{\mathrm{R}}=L_t^{\mathrm{R}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}\text{(C7)}& {\displaystyle }{L}_t^{\mathrm{C}}=\sum _{v=\mathrm{1}}^{v_{max}}{L}_{v,t}^{\mathrm{C}},\text{(C8)}& {\displaystyle }{L}_{v+\mathrm{1},t+\mathrm{1}}^{\mathrm{C}}=L_{v,t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v,t}^{\mathrm{C},\mathrm{H}},v<v_{max}-\mathrm{1}\text{(C9)}& {\displaystyle }\begin{array}{rl}{L}_{v_{max},t+\mathrm{1}}^{\mathrm{C}}& =L_{v_{max},t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v_{max},t}^{\mathrm{C},\mathrm{H}}-{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}}\\ & +L_{v_{max}-\mathrm{1},t}^{\mathrm{C}}-{\mathrm{\Delta }}_{v_{max}-\mathrm{1},t}^{\mathrm{C},\mathrm{H}},\end{array}\text{(C10)}& {\displaystyle }{L}_{\mathrm{1},t+\mathrm{1}}^{\mathrm{C}}={\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}}\end{array}$$

$$\begin{array}{l}{\displaystyle }{\mathrm{\Delta }}_{v,t}^{\mathrm{C},\mathrm{H}}\le L_{v,t}^{\mathrm{C}},v<v_{max}\\ {\displaystyle }{\mathrm{\Delta }}_{v_{max},t}^{\mathrm{C},\mathrm{H}}+{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}}\le L_{v_{max},t}^{\mathrm{C}}\\ {\displaystyle }{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}}=\sum _{v=\mathrm{1}}^{v_{max}}{\mathrm{\Delta }}_{v,t}^{\mathrm{C},\mathrm{H}}-{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}}\end{array}$$

C1.2 Fossil fuels

$$\begin{array}{}\text{(C11)}& X_{t+\mathrm{1}}^{\mathrm{F}}=X_t^{\mathrm{F}}-{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}-{\mathrm{\Delta }}_t^{\mathrm{F},n}+{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{D}}\end{array}$$

C1.3 Other primary resources

$$\begin{array}{}\text{(C12)}& X_t^{\mathrm{O}}=X_{\mathrm{0}}^{\mathrm{O}}\left[{\mathit{\alpha }}^{\mathrm{o},\mathrm{l}}{\displaystyle \frac{{\mathrm{\Pi }}_t}{{\mathrm{\Pi }}_{\mathrm{0}}}}+\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{o},\mathrm{l}}\right)(\mathrm{1}+{\mathit{\kappa }}^{\mathrm{o},\mathrm{2}}{)}^t\right]\end{array}$$

C1.4 Intermediate products

$$\begin{array}{}\text{(C13)}& {\displaystyle }{x}_t^p={\mathit{\theta }}_t^p{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}\text{(C14)}& {\displaystyle }{x}_t^n={\mathit{\theta }}^n{\mathrm{\Delta }}_t^{\mathrm{F},n}\text{(C15)}& {\displaystyle }{x}_t^n=x_t^{n,\mathrm{c}}+x_t^{n,b\mathrm{2}}\text{(C16)}& {\displaystyle }{x}_t^{\mathrm{c}}={\mathit{\theta }}_t^{\mathrm{c}}{\left({\mathit{\alpha }}^n{\left(L_t^{\mathrm{A},\mathrm{c}}\right)}^{{\mathit{\rho }}_n}+\left(\mathrm{1}-{\mathit{\alpha }}^n\right){\left(x_t^{n,\mathrm{c}}\right)}^{{\mathit{\rho }}_n}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_n}}\text{(C17)}& {\displaystyle }\begin{array}{rl}{x}_t^{\mathrm{c},\mathrm{b}\mathrm{2}}& ={\mathit{\theta }}_t^{\mathrm{c},\mathrm{b}\mathrm{2}}\left({\mathit{\alpha }}^n{\left(L_t^{\mathrm{A},\mathrm{b}\mathrm{2}}\right)}^{{\mathit{\rho }}_n}\right.\\ & {\left.+\left(\mathrm{1}-{\mathit{\alpha }}^n\right){\left(x_t^{n,\mathrm{b}\mathrm{2}}\right)}^{{\mathit{\rho }}_n}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_n}}\end{array}\text{(C18)}& {\displaystyle }{x}_t^{\mathrm{b}\mathrm{1}}={\mathit{\theta }}^{\mathrm{b}\mathrm{1}}{x}_t^{\mathrm{c},\mathrm{b}}\text{(C19)}& {\displaystyle }\begin{array}{rl}{x}_t^{\mathrm{b}\mathrm{2}}& ={\mathit{\theta }}^{\mathrm{b}\mathrm{2}}\\ & {\left({\left({\mathit{\alpha }}^{\mathrm{b}\mathrm{2}}\right)}^{{\mathit{\theta }}_t^{\mathrm{b}\mathrm{2},K}}{\left(K\right)}^{{\mathit{\rho }}_{b\mathrm{2}}}+\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{b}\mathrm{2}}\right){\left(x_t^{\mathrm{c},\mathrm{b}\mathrm{2}}\right)}^{{\mathit{\rho }}_{\mathrm{b}\mathrm{2}}}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_{b\mathrm{2}}}}\end{array}\text{(C20)}& {\displaystyle }{x}_t^{\mathrm{l}}={\mathit{\theta }}^{\mathrm{P}}{\left({\mathit{\alpha }}^{\mathrm{l}}{\left(L_t^{\mathrm{P}}\right)}^{{\mathit{\rho }}_{\mathrm{l}}}+\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{l}}\right){\left(x_t^{\mathrm{c},\mathrm{l}}\right)}^{{\mathit{\rho }}_{\mathrm{l}}}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{l}}}}\text{(C21)}& {\displaystyle }{x}_t^{\mathrm{w}}=\sum _{v=\mathrm{1}}^{v_{max}}{\mathit{\theta }}_{v,t}^{\mathrm{w}}\mathrm{\Delta }{}_{v,t}^{\mathrm{C},\mathrm{H}}\end{array}$$

C1.5 Final goods and services

$$\begin{array}{}\text{(C22)}& {\displaystyle }{Y}_t^{\mathrm{e}}={\mathit{\theta }}_t^{\mathrm{e}}{\left({\mathit{\alpha }}^{\mathrm{e}}{\left(x_t^{\mathrm{b}\mathrm{1}}\right)}^{{\mathit{\rho }}_e}+\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{e}}\right){\left(x_t^p+x_t^{\mathrm{b}\mathrm{2}}\right)}^{{\mathit{\rho }}_e}\right)}^{\frac{\mathrm{1}}{{\mathit{\rho }}_e}}\text{(C23)}& {\displaystyle }{Y}_t^{\mathrm{l}}={\mathit{\theta }}_t^{\mathrm{l}}{x}_t^{\mathrm{l}},\text{(C24)}& {\displaystyle }{Y}_t^{\mathrm{w}}={\mathit{\theta }}_t^{y_{\mathrm{w}}}{x}_t^{\mathrm{w}}\text{(C25)}& {\displaystyle }\begin{array}{rl}{Y}_t^{\mathrm{r}}& ={\mathit{\theta }}^{\mathrm{r}}\left[\sum _{i=\mathrm{A},\mathrm{P},\mathrm{C}}{\mathit{\alpha }}^{i,\mathrm{r}}{\left(L_t^i\right)}^{{\mathit{\rho }}_r}\right.\\ & {\left.+\left(\mathrm{1}-\sum _{i=\mathrm{A},\mathrm{P},\mathrm{C}}{\mathit{\alpha }}^{i,r}\right){\left(L_t^{\mathrm{N}}+{\mathit{\theta }}^{\mathrm{R}}{L}_t^{\mathrm{R}}\right)}^{{\mathit{\rho }}_r}\right]}^{\frac{\mathrm{1}}{{\mathit{\rho }}_r}}\end{array}\text{(C26)}& {\displaystyle }\begin{array}{rl}& Y_t^{\mathrm{o}}={\mathit{\theta }}_t^{o,\mathrm{1}}\\ & \left[\begin{array}{c}{X}_t^{\mathrm{O}}-\frac{\mathrm{1}}{{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{o}}}[c^{\mathrm{o},\mathrm{c}}\frac{x_t^{\mathrm{c}}}{A_t}+c^{\mathrm{o},\mathrm{cb}}\frac{x_t^{\mathrm{c},\mathrm{b}\mathrm{2}}}{{\mathit{\theta }}_t^{\mathrm{c},\mathrm{b}\mathrm{2}}}+c^{\mathrm{o},\mathrm{f}}\frac{Y_t^{\mathrm{f}}}{{\mathit{\theta }}_t^{\mathrm{f}}}+c^{\mathrm{o},\mathrm{p}}{x}_t^{\mathrm{p}}+c^{\mathrm{o},n}{x}_t^n+c^{\mathrm{o},\mathrm{b}}{x}_t^{\mathrm{b}\mathrm{1}}\\ +c^{\mathrm{o},\mathrm{b}\mathrm{2}}{x}_t^{\mathrm{b}\mathrm{2}}+c^{\mathrm{o},\mathrm{l}}{x}_t^{\mathrm{l}}+c^{\mathrm{o},y_{\mathrm{l}}}\frac{{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{l}}{Y}_t^{\mathrm{l}}}{{\mathit{\theta }}_t^{\mathrm{l}}}+c_t^{\mathrm{o},\mathrm{w}}{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{H}}+c^{\mathrm{o},y_{\mathrm{w}}}{x}_t^{\mathrm{w}}\\ +c^{\mathrm{o},\mathrm{r}}{L}_t^{\mathrm{R}}+c^{\mathrm{p}}{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{C}}+C_t^{\mathrm{N}}+C_t^{\mathrm{R}}+C_t^{\mathrm{F}}+C_t^{\mathrm{H}}+C_t^{\mathrm{P}}+C_t^{\mathrm{C},\mathrm{N}}]\end{array}\right]\end{array}\text{(C27)}& {\displaystyle }\begin{array}{rl}{\mathit{y}}_t& =\left(y_t^{\mathrm{f}},y_t^{\mathrm{l}},y_t^{\mathrm{e}},y_t^{\mathrm{w}},y_t^{\mathrm{r}},y_t^{\mathrm{o}}\right)\\ & =\left(Y_t^{\mathrm{f}},Y_t^{\mathrm{l}},Y_t^{\mathrm{e}},Y_t^{\mathrm{w}},Y_t^{\mathrm{r}},Y_t^{\mathrm{o}}\right)/{\mathrm{\Pi }}_t\end{array}\end{array}$$

C1.6 Technology (deterministic)

$$\begin{array}{}\text{(C28)}& {\displaystyle }{A}_t={\displaystyle \frac{A_T{A}_{\mathrm{0}}{e}^{{\mathit{\kappa }}_{\mathrm{c}}t}}{A_T+A_{\mathrm{0}}\left(e^{{\mathit{\kappa }}_{\mathrm{c}}t}-\mathrm{1}\right)}}\text{(C29)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{v,t}^{\mathrm{w}}=\\ & \left\{\begin{array}{l}\mathrm{0.00001}\text{ if }v\le \underset{\mathrm{¯}}{v}\\ {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_v^{\mathrm{w}}(\mathrm{1}+{\mathit{\kappa }}_v^{\mathrm{w}}t)\text{ if }v>\underset{\mathrm{¯}}{v}\end{array}\right.,{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_v^{\mathrm{w}}=\exp \left({\mathit{\psi }}_a-{\displaystyle \frac{{\mathit{\psi }}_b}{(v-\underset{\mathrm{¯}}{v})}}\right)\end{array}\text{(C30)}& {\displaystyle }{\mathit{\theta }}_t^i={\mathit{\theta }}_{\mathrm{0}}^i(\mathrm{1}+{\mathit{\kappa }}^i{)}^t,i=f,e,l,y^{\mathrm{w}},o\end{array}$$

C1.7 Technology (stochastic)

$$\begin{array}{}\text{(C31)}& A_t={\displaystyle \frac{A_T(J_{\mathrm{1},t},J_{\mathrm{2},t})A_{\mathrm{0}}{e}^{{\mathit{\kappa }}_{\mathrm{c}}t}}{A_T(J_{\mathrm{1},t},J_{\mathrm{2},t})+A_{\mathrm{0}}\left(e^{{\mathit{\kappa }}_{\mathrm{c}}t}-\mathrm{1}\right)}}\end{array}$$

C1.8 Costs

$$\begin{array}{}\text{(C32)}& {\displaystyle }{C}_t^{\mathrm{N},\mathrm{A},\mathrm{R}}={\mathit{\xi }}_{\mathrm{0}}^n\left({\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}\right)+{\mathit{\xi }}_{\mathrm{1}}^n{\left({\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}}+{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}\right)}^{\mathrm{2}}\text{(C33)}& {\displaystyle }{C}_t^{N,R}={\mathit{\xi }}_{\mathrm{0}}^{\mathrm{R}}{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}+{\mathit{\xi }}_{\mathrm{1}}^{\mathrm{R}}{\left({\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}}\right)}^{\mathrm{2}}\text{(C34)}& {\displaystyle }{C}_t^F={\mathit{\xi }}_{\mathrm{1}}^F{\left({\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}}+{\mathrm{\Delta }}_t^{F,n}\right)}^{\mathrm{2}}\left({\displaystyle \frac{X_{\mathrm{0}}^F+{\mathrm{\Delta }}^{\mathrm{F},\mathrm{D}}}{X_t^F+{\mathrm{\Delta }}^{F,D}}}\right)\text{(C35)}& {\displaystyle }{C}_t^H={\mathit{\xi }}_{\mathrm{0}}^H{\left({\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{H}}-{\mathrm{\Delta }}_t^{C,C}\right)}^{\mathrm{2}}+\sum _v{\displaystyle \frac{{\mathit{\xi }}_{\mathrm{1}}^H}{L_{v,t+\mathrm{1}}^{\mathrm{C}}+{\mathit{\xi }}_{\mathrm{2}}^H}}\text{(C36)}& {\displaystyle }{C}_t^{A,P}={\mathit{\xi }}_{\mathrm{1}}^{\mathrm{P}}{\left({\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}}\right)}^{\mathrm{2}}\text{(C37)}& {\displaystyle }{C}_t^{C,N}={\mathit{\xi }}_{\mathrm{0}}^{C,N}{\mathrm{\Delta }}_t^{C,N}+{\mathit{\xi }}_{\mathrm{1}}^{C,N}{\left({\mathrm{\Delta }}_t^{C,N}\right)}^{\mathrm{2}}\end{array}$$

C1.9 Preferences

$$\begin{array}{}\text{(C38)}& {\displaystyle }u\left(\mathit{y}\right)={\displaystyle \frac{\mathcal{C}{\left(\mathit{y}\right)}^{\mathrm{1}-\mathit{\gamma }}}{\mathrm{1}-\mathit{\gamma }}}\text{(C39)}& {\displaystyle }\begin{array}{rl}& \log \left(\mathcal{C}\left(\mathit{y}\right)\right)=\\ & \sum _{q=f,l,e,w,r,o}\left({\displaystyle \frac{{\mathit{\alpha }}_q+{\mathit{\beta }}_q\mathcal{C}\left(\mathit{y}\right)}{\mathrm{1}+\mathcal{C}\left(\mathit{y}\right)}}\right)\log \left(y_t^q-{\underset{\mathrm{¯}}{y}}^q\right)\end{array}\end{array}$$

C1.10 Population

$$\begin{array}{}\text{(C40)}& {\mathrm{\Pi }}_t={\displaystyle \frac{{\mathrm{\Pi }}_T{\mathrm{\Pi }}_{\mathrm{0}}{e}^{{\mathit{\kappa }}^{\mathit{\pi }}t}}{{\mathrm{\Pi }}_T+{\mathrm{\Pi }}_{\mathrm{0}}\left(e^{{\mathit{\kappa }}^{\mathit{\pi }}t}-\mathrm{1}\right)}}\end{array}$$

C1.11 Welfare

$$\begin{array}{}\text{(C41)}& \mathrm{\Omega }=\mathbb{E}\left\{\sum _{t=\mathrm{0}}^{\mathrm{\infty }}{\mathit{\delta }}^t\mathcal{U}({\mathit{S}}_t,{\mathit{a}}_t)\right\}.\end{array}$$

with $\mathcal{U}({\mathit{S}}_t,{\mathit{a}}_t)=u\left({\mathit{y}}_t\right){\mathrm{\Pi }}_t$, $\mathit{S}:=(\mathbf{L},\mathbf{X},\mathbf{J})$, and

$$\begin{array}{rl}{\mathit{a}}_t& =({\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{A}},{\mathrm{\Delta }}_t^{\mathrm{N},\mathrm{R}},{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{N}},{\mathrm{\Delta }}_t^{\mathrm{A},\mathrm{P}},{\mathrm{\Delta }}_t^{\mathrm{C},\mathrm{A}},{\mathrm{\Delta }}_{\mathrm{1},t}^{\mathrm{C},\mathrm{H}},\\ & \mathrm{\dots },{\mathrm{\Delta }}_{v_{max},t}^{\mathrm{C},\mathrm{H}},{\mathrm{\Delta }}_t^{C,C},{\mathrm{\Delta }}_t^{\mathrm{F},\mathrm{E}},{\mathrm{\Delta }}_t^{\mathrm{F},n},L_t^{A,F},L_t^{\mathrm{A},\mathrm{B}},{\mathit{x}}_t,{\mathit{y}}_t).\end{array}$$

Figure C1 Changes in potential crop yields under Representative Concentration Pathway (RCP) 6 scenario in 2100.

Figure C2 Consumption of fertilizers and biofuels. The solid lines represent the deterministic–baseline solutions. Squares represent the deterministic–optimistic solutions or their distances to the deterministic-baseline solutions. Marks represent the deterministic–pessimistic solutions or their distances to the deterministic-baseline solutions. The shaded areas are the ranges of the sample paths' distances to the deterministic-baseline solutions. Pluses represent the average paths, the dashed lines represent the 10 % quantile, circles represent the median paths, and the dashed–dotted lines represent the 90 % quantile. Upward-pointing triangles and downward-pointing triangles represent the distances of the stochastic–optimistic solutions and the stochastic–pessimistic solutions, respectively, to the deterministic-baseline solutions.

Table C1 Model exogenous variables.

Table C2 Model endogenous variables.

Table C3 Baseline parameters.