Understanding Uncertainty in Market-Mediated Responses to US Oilseed Biodiesel Demand: Sensitivity of ILUC Emission Estimates to GLOBIOM Parametric Uncertainty

Understanding Uncertainty in Market-Mediated Responses to US Oilseed Biodiesel Demand: Sensitivity of ILUC Emission Estimates to GLOBIOM Parametric Uncertainty Neus Escobar,*Hugo Valin, Stefan Frank,*Diana Galperin, Christopher M. Wade, Leopold Ringwald, Daniel Tanner, Niklas Hinkel, Petr Havlík, Justin S. Baker, Sharyn Lie, and Christopher Ramig Cite This: https://doi.org/10.1021/acs.est.3c09944 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: The life cycle greenhouse gas (GHG) emissions of biofuels depend on uncertain estimates of induced land use change (ILUC) and subsequent emissions from carbon stock changes. Demand for oilseed-based biofuels is associated with particularly complex market and supply chain dynamics, which must be considered. Using the global partial equilibrium model GLOBIOM, this study explores the uncertainty in market-mediated impacts and ILUC-related emissions from increasing demand for soybean biodiesel in the United States in the period 2020−2050. A one-ata-time (OAT) analysis and a Monte Carlo (MC) analysis are performed to assess the sensitivity of modeled ILUC-GHG emissions intensities (gCO2e/MJ) to varying key economic and biophysical model parameters. Additionally, the influence of the approach on the simulation of future ILUC effects is explored using two alternative ILUC-GHG metrics: a comparative-static approach for 2030 and a recursive-dynamic approach using model outputs through 2050. We find that projected ILUC-GHG values largely vary based on which vegetable oils replace diverted soybean oil, market responses to coproducts, and the carbon content of land converted for agricultural use. These are all, in turn, subject to decision uncertainty through the choice of the modeling approach and the time horizon considered for each ILUC-GHG metric. Given the longer simulation period, ILUC-GHG emission uncertainty ranges increase under the recursive-dynamic approach (42.4 ±25.9 gCO2e/MJ) compared to the comparative-static approach (40.8 ±20.5 gCO2e/MJ). The combination of MC analysis with other techniques such as Bayesian Additive Regression Trees (BART) is powerful for understanding model behavior and clarifying the sensitivity of market responses, ILUC, and associated GHG emissions to specific model parameters when simulated with global economic models. The BART reveals that biophysical parameters generate more linear ILUC-GHG responses to changes in assumed parameter values while changes in economic parameters lead to more nonlinear ILUC-GHG results as multiple effects at the interplay of food, feed, and fuel uses overlap. The choice of the recursive-dynamic metric allows capturing the longer-term evolution of ILUC while generating additional uncertainties derived from the baseline definition. KEYWORDS: biofuels, climate change, economic model, greenhouse gas, international trade, life cycle assessment, spillover, transport 1. INTRODUCTION The potential contribution of biofuels to climate change mitigation has received significant attention in the literature. 1−3 The estimation of land use change (LUC) emissions associated with biofuels has been debated intensively, motivated partially by the fact that potential future LUC impacts cannot be measured but only modeled. 4−8 LUC can refer to either land conversion to grow biofuel feedstocks (often referred to as direct LUC) or the subsequent land transformation to support other uses globally, such as food and feed (referred to as indirect LUC). Induced land use change (ILUC) hereinafter designates the combination of the two components, capturing the overall market-mediated displacement of land uses in response to an increased demand for biomass for fuel. 9,10 ILUC results in net greenhouse gas (GHG) emissions when global net land carbon stocks are lost through successive land conversions. 11,12 The two methods commonly used to quantify the GHG emissions intensity of biofuels are attributional life cycle assessment (ALCA) and consequential LCA (CLCA). ALCA considers impacts from feedstock production up to combustion, often including direct LUC. 13−16 Quantifying indirect LUC emissions requires economic reasoning under conReceived: May 29, 2024 Revised: November 27, 2024 Accepted: December 2, 2024 Published: December 18, 2024 Articlepubs.acs.org/est © 2024 The Authors. Published by American Chemical Society A https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX This article is licensed under CC-BY-NC-ND 4.0 Downloaded via 62.42.100.239 on December 20, 2024 at 20:00:16 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. sequential approaches. 17,18 For biofuels, CLCA aims to represent how increased demand for crops causes agricultural land expansion, crop displacement, intensification, and changes to crop and livestock production emissions, depending on relative productivity and price adjustments across sectors and regions. 19,20 CLCA often applies economic modeling to determine the extent and location of land conversion and the ecosystems and land uses affected. 21−23 In a review of the science of biofuel LCA, the US National Academies of Science, Engineering, and Medicine found that CLCA is preferable when seeking to understand the consequences of decisions or actions that alter overall quantities of biofuel consumed. 24 Recent ALCA studies show that vegetation loss remains a large contributor to the GHG emission intensity of biofuel feedstocks, showing that direct LUC emissions can be higher than ILUC emission values for biofuels produced in specific locations when using spatially explicit carbon stock and land use data. 25−27 Global economic models provide consistent frameworks to assess ILUC impacts of biofuels under CLCA approaches. 28−30 These models capture key economic mechanisms that shape the global distribution of land uses, such as yield and cropland area responses to changes in land availability and trade. 31 For instance, the Computable General Equilibrium (CGE) model GTAP-BIO 32,33 includes several food and nonfood biofuel sectors −and their coproducts−to simulate economy-wide effects and land cover change due to expanded biofuel production. GLOBIOM 34,35 is a Partial Equilibrium (PE) model representing agriculture and forestry sectors, with a spatially explicit land use representation. Both models have been widely used for the analysis of biofuel policies. 32,36−40 ILUC impacts of biofuels have been analyzed extensively, initially for the so-called first-generation biofuels 41−43 and more recently for nonfood feedstocks. 44−47 Particularly, the ILUC emission estimation for feedstocks for aviation biofuel production as part of the Carbon Offsetting and Reduction Scheme for International Aviation (CORSIA) uses GTAP-BIO and GLOBIOM to propose a harmonized ILUC emission intensity value per feedstock, including nonfood crops and residues. 48 In the context of the United States (US), studies have mainly focused on corn ethanol 42,43,49 and soybean biodiesel, 5,50,51 with the impacts of oilseed production and trade recently gaining more attention. 52−55 As international demand for oilseeds and vegetable oil for various uses continues to increase worldwide, there are concerns that these increasing demands impact forest and natural vegetation loss, especially in the tropics. 56−59 Global combined biodiesel and renewable diesel production expanded from around 20.2 billion liters in 2010 to around 52.7 billion liters in 2021. 60 Following a 5-fold increase in US biodiesel supply between 2010 to 2021, the US is now the second largest biodiesel producer behind the EU, accounting for 20% of global production. 60−63 Most US biodistillate fuel production over this period has been sourced from vegetable oils, with soybean oil representing the plurality of biodistillate feedstock. Estimating ILUC emission intensities of vegetable oil-based fuels has diverse sources of uncertainty including the choice of the modeling framework and accompanying assumptions to quantify current and future effects, e.g., baseline year, analytic horizon (decision uncertainty) and the lack of full understanding of the underlying complex land use dynamics (epistemic uncertainty). 42 The latter is related to the input parameters and specific values needed to simulate the processes being modeled and is also referred as parametric uncertainty. 64,65 CLCA modeling requires numerous data and input parameters, which represent and influence both economic and biophysical dynamics, to capture responses across interlinked food, feed, and fuels markets. Key dynamics for evaluating ILUC effects include cropland intensification/ extensification as well as impacts on land use and management and livestock production. 66 Previous work found that, among the parameters determining long-term ILUC effects, the most decisive relate to costs of land transformation, endogenous productivity responses, and cropland availability. 66 Past uncertainty analyses of biofuels focused on the effects of alternative yield elasticities to crop prices, demand elasticities, trade elasticities, and land transformation (expansion) elasticities to land rents. 37,40,67 Studies performing sensitivity analyses with the CGE model GTAP-BIO 32 highlight the role of the elasticity governing land conversion between cropland, pastureland, and managed forestland. 37,68 When assessing parametric uncertainty in GTAP-BIO with MC analysis, 64 the yield elasticity to price, the Armington trade elasticities, the GHG emissions factor (EF) for cropland-to-pasture conversion, and the productivity of newly converted cropland were found to contribute most of the variance in ILUC emission values for US soybean biodiesel. The size and complexity of global economic models pose challenges for designing, implementing, running, and interpreting large-scale sensitivity analyses intended to propagate uncertainty from parameters to model outcomes. As such, only a few studies have devoted efforts to execute these methods. 37,44,49,64,65 Using the GLOBIOM global economic modeling framework, 69 this study aims to assess the sensitivity of two different measures of ILUC-GHG emissions intensities (hereinafter referred to as “ILUC-GHG values”) of US soybean biodiesel to the uncertainty in key model parameters influencing ILUC responses. This is done by combining a oneat-a-time (OAT) analysis, where each of the selected parameters is varied individually over a predetermined range; and a Monte Carlo (MC) analysis, where all parameters are given specific probability distributions and varied randomly and simultaneously. Additionally, we assess the influence of the choice of the time horizon and ILUC emission accounting procedure by using two different approaches to derive future ILUC-GHG values: a comparative-static approach, in which modeled changes in 2030 are amortized over 25 years of biodiesel production, and a recursive-dynamic approach using model outputs through 2050. Thus, the goal of the study is to better understand GLOBIOM model behavior and the role of certain parameters and assumptions when estimating ILUCGHG values. Overall, this study contributes to the understanding of the factors driving market-mediated responses to increased demand for biobased products, while identifying areas of future research for a better representation of oilseed biofuels in global economic models. 2. METHODS 2.1. Modeling Framework. The global recursive-dynamic PE model GLOBIOM 40,69 is applied to assess the effects of an increased demand for soybean biodiesel in the US in the period 2020−2050, hereinafter called shock. GLOBIOM computes a global equilibrium in agricultural and forest product markets in 10-year time steps through the period 2000−2050 by choosing land use and processing activities that maximize welfare subject to resource, technological, demand, Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX B and policy constraints. The model does not represent marketmediated effects of biofuel policies on sectors of the economy other than land use (forestry and agriculture, including crop and livestock). For this analysis, the model is aggregated into 57 world regions. 34,36 The baseline represents economic developments according to the Shared Socioeconomic Pathway 2 scenario, 70 while historical yields were recalibrated so that they include 2017−2019 data. 71 Vegetable oil prices were also recalibrated to align with relative price projections both for the US versus the rest of the world (ROW), respectively based on data from the United States Department of Agriculture (USDA) and the Food and Agriculture Organization (FAO) Further details are included in Section S1.1 of the electronic Supporting Information (ESM). 2.2. Scenarios and ILUC Metrics. This study models a soybean biodiesel consumption shock in the US introduced in 2021 to progressively reach a total additional demand of 126.9 PJ/year in 2030, which corresponds to 1 billion gallons of gasoline equivalent, or roughly a doubling of current US consumption. This is compared to a baseline that keeps global biofuel volumes constant through 2050 at 2020 levels. To understand the uncertainty brought about by the ILUC emission accounting procedure, two alternative methods are applied to calculate ILUC-GHG values, hereinafter referred to as comparative-static and recursive-dynamic (see equations in the Appendix). The former projects the ILUC-GHG values in 2030 −the year in which the full biofuel mandate is reached − over a 10-year period. This approach tries to isolate the effects of the biofuel shock in the short term, assuming the impacts will be equally distributed throughout the amortization period of 25 years. Alternatively, the recursive-dynamic approach considers a constant demand of 126.9 PJ/year above baseline levels until 2050, thereby capturing market and land use developments in the longer term. This corresponds to biofuel volumes equivalent to 25 years of total mandated biodiesel consumption in 2020−2050 (Figure S1 in ESM), in line with the 25-year amortization period under the comparative-static approach. Under both approaches, ILUC-GHG values are calculated considering CO2e emissions/removals from changes across carbon pools according to IPCC guidelines (IPCC 2006) (eq 1). This refers to net changes in aboveand below-ground biomass, dead wood, litter, and harvested wood products, while additional emissions arise from peatland oxidation (eq 2) and changes in soil organic carbon (SOC) relative to reference soils. In GLOBIOM, net changes in carbon stocks are the result of the land transitions simulated, ultimately leading to natural vegetation loss, natural vegetation reversion, and cropland expansion. The comparative-static ILUC-GHG value considers changes in land use-related emissions for 2030 when the biofuel consumption mandate is reached (eq 3), whereas the recursive-dynamic ILUC-GHG value considers cumulative emissions up to 2050 (eq 4). The comparativestatic ILUC-GHG value is often used in scientific literature looking at short-term responses of a biofuel shock 32,39,72 and hence facilitates comparison of our results, while the recursivedynamic ILUC-GHG value gives additional insights on how the biofuel emissions impacts evolve over time. 49,73 2.3. Sensitivity Analysis. Sensitivity analysis is applied to identify how variation in input parameter values causes variation in output variables. In this case, sensitivity in ILUC-GHG values is assessed by changing input parameters incrementally and individually, as well as stochastically and simultaneously, in OAT and MC analyses, respectively. The two analyses are combined to understand the model’s performance for ILUC-GHG value estimation, shedding light on the role of the varied parameters in determining uncertainty. Eleven model parameters and their associated probability distributions were chosen based on expert knowledge and experience from previous GLOBIOM studies. 40,74,75 They include seven economic and four biophysical parameters (Table S1 in ESM). To initially understand the influence of these parameters on estimated ILUC-GHG values, the OAT analysis varies each parameter alone while keeping the rest at the central value, i.e., at the model default. For each key parameter, eight alternative values were considered, four below and four above the central value, covering the parameter ranges detailed in ESM Table S1. For each of the 99 different parameter combinations (eight alternatives for all 11 parameters, and a central case with all default values), baseline and shock scenarios were run to calculate both ILUC-GHG metrics considered. MC simulation is a stochastic technique that produces a range of results from which a probability distribution of modeled results can be inferred. MC analysis involves solving the model across a number of runs that is sufficiently high relative to the sample size used for input parameters. 76 In each run, values for specific model input parameters are randomly selected from defined probability distributions for each stochastic input parameter. In this study, 1,000 runs are performed for both the shock and baseline scenarios to estimate the difference in model outcomes for each scenario pair with the same combination of input parameters. In contrast with the OAT analysis, several parameter values are drawn independently across regions (substitution elasticity− vegetable oils), commodities (trade elasticity−vegetable oils), or land types (land expansion into natural vegetation) for each MC run. This yields 72 individual probability distributions in total, considering the 11 parameters varied and the specific regions and products to which they apply. Lastly, to identify the role of each parameter in MC analysis, Bayesian Additive Regression Trees (BART) analysis 77 was applied ex-post to approximate the functional form of GLOBIOM outcomes with respect to the parameters under investigation. More specifically, the BART model is trained on the sampled economic and biophysical shock values as covariates, in order to explain the variation in ILUC-GHG values obtained across all MC runs. More details on the approach are detailed in the ESM (Section S1.6). 3. RESULTS 3.1. Central Case. This section presents the central case (with default model parameter values) to contextualize the subsequent sensitivity analysis and summarize important model dynamics in response to the shock. GLOBIOM estimates that the US biodiesel consumption shock leads to three key market-mediated effects that determine the ILUC effects and related GHG intensity values in the central case: 1) global vegetable oil demand for nonfuel uses decreases, 2) other vegetable oils (mainly palm oil) substitute for soybean oil for nonfuel uses globally, and 3) the regional distribution of soybean production shifts at the margin. The first and second effects are integral to understanding the ILUC implications of a US soybean biodiesel shock in the central case. The US sources most of the additional soybean oil needed to supply the biofuel shock (a total of 3.5 Mt) in 2030 through Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX C increased domestic production (1.9 Mt), reduction of net trade (−1.2 Mt), and domestic reduction of nonfuel uses of soybean oil (−0.3 Mt) (Figure 1). This significantly affects global vegetable oil markets and land use dynamics in non-US regions (Figure S2). The increased US soybean oil demand and the decreased exports raise global prices for soybean oil and, subsequently, other vegetable oils (by 6% and 3% in 2030, respectively). While global supply of vegetable oils increases in 2030 (2.1 Mt), demand for food and other nonfuel applications decreases (−1.3 Mt); both because vegetable oil prices rise, especially for food uses, and because vegetable oils are not perfect substitutes. Palm and rapeseed oils partly compensate for decreases in soybean oil consumption (−3.1 Mt) for nonfuel applications globally, with additional 1.8 Mt consumed in 2030. This strongly affects GHG emissions from ILUC in GLOBIOM, as additional palm oil originates from Southeast Asia (SEA), where palm expansion in the past entailed rainforest conversion and peatland oxidation. As the US increases domestic soybean production in response to the shock (+5.7 Mt, + 1.4 Mha harvested area in 2030), domestic crushing increases availability of soybean meal (+7.7 Mt), demanded for animal feeding in international markets. While prices of soybean oil increase 27.1% in the US and 5.5% globally in 2030, prices of soybean meal decrease compared to the baseline by 7.8% domestically and 3.3% globally. The increased availability of relatively cheap US soybean meal for livestock feeding puts pressure on less productive producers in Asia but also reduces demand for soybeans and soybean meal from other leading exporters, mainly Argentina and Brazil, who, together, account for around 60% of the global market and exports of soybeans. As world average soybean prices decrease, soybean area expansion in South America (SAM) in 2020−2030 is more limited and declines by 1.4 percentage points, to 33.8%, compared to 35.2% expansion in the baseline. This translates into lower cropland expansion in SAM (+0.3 Mha, + 0.28% in 2030) with the shock compared to the baseline. In 2050, cropland expands in SAM by 0.1 Mha (+30.8% in 2020−2050 vs +30.7% in the baseline) due to increased demand for other agri-food crops globally, mainly through the livestock rebound effect, which herein describes the increased profitability of livestock sectors globally given the increased availability of relatively cheap meals in the market, and the associated increase in animal products consumption. As a result, animal production systems intensify (especially pig, poultry, and dairy), demanding more feed concentrates and other grains to complement feed rations. Global cereal demand for animal feed increases by 1.3 Mt as well as related harvested areas (+0.2 Mha in 2030) (Figure 1). Grassland areas also expand in regions with relatively extensive livestock production systems, e.g., US, SAM, and Southern Asia (SAS). ILUC effects increase global GHG emissions relative to the baseline. The comparative-static ILUC-GHG value is estimated at 29.4 gCO2e/MJ (Figure 1) with the most important emissions source being peatland oxidation (36.0 gCO2e/MJ), followed by natural land conversion (16.3 gCO2e/MJ), and forgone natural land reversion (2.4 gCO2e/MJ). These emissions are partly compensated by enhanced carbon sequestration from agricultural biomass (−22.6 gCO2e/MJ) and SOC (−2.7 gCO2e/MJ). Palm expansion in SEA contributes significantly to ILUC-related emissions, especially through natural land conversion and peatland oxidation. Simultaneously, palm plantations sequester more carbon per hectare than annual crops thereby generating some carbon sink in agricultural biomass. The size of the total biomass effect from palm depends on the type of land that palm plantations are expanding onto. If palm plantations replace rainforests, there is net loss in total biomass carbon. Since the model estimates greater expansion onto previously undisturbed peatlands, peatland emissions associated with palm expansion outweigh the biomass sequestration effect. These results in the central case scenario are higher than the GTAP-BIO value of around 20 gCO2e/MJ for US soybean biodiesel. 32,55,78 Notable differences in assumptions between these studies include biofuel volumes, shock year, and biofuel processing efficiencies. However, as discussed further below, this 20 gCO2e/MJ result falls within the range of GLOBIOM results described in our sensitivity analysis. Figure 1. Absolute changes in global harvested area (left) and global land cover (center) when comparing the shock scenario to the baseline scenario in 2030 and 2050; and corresponding ILUC-GHG values under comparative-static and recursive-dynamic approaches (right), respectively, broken down by GHG source. Each “x” represents the net value across categories depicted. Different multicropping intensities across crops and world regions explain the difference in total harvested area vs total cropland changes. Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX D Global GHG emissions increase over time, primarily due to the increasing production and use of palm oil to substitute for decreased US soybean oil availability for food and other nonfuel uses (Figure S2 in ESM). The recursive-dynamic ILUC-GHG value, which, in contrast to the comparative-static value, accounts for these lagged effects of the shock as the model recursively solves timesteps through 2050, is 33.6 gCO2e/MJ (Figure 1). Most emissions come from peatland oxidation (38.8 gCO2e/MJ) while natural land conversion and associated biomass losses contribute with 9.6 gCO2e/MJ. Agricultural biomass growth offsets emissions by 25.4 gCO2e/ MJ, mostly through increases in oil palm areas followed by corn and other cereals for the livestock rebound effect (Figure 1) (see eq 1 in Appendix). 3.2. Sensitivity Analysis. 3.2.1. One-at-a-Time Sensitivity Analysis. In the OAT analysis, some, but not all, of the parameter alterations produce considerable variation in the ILUC-GHG values relative to the central case (Figure 2). For the comparative-static ILUC-GHG value (mean ±standard deviation, gCO2e/MJ), the largest variations are associated with the trade elasticity of vegetable oils (35.8 ±13.6), the expansion response of palm into peatland (29.4 ±12.3), and the EF for carbon sequestration in biomass in palm plantations (29.9 ±12.3). These results underpin the findings of the central case and highlight the importance of land use spillovers in the SEA and SAM regions as major contributors to ILUC-GHG values and their uncertainty. The demand elasticity for vegetable oils (38.3 ±10.7) and the substitution elasticity among vegetable oils (33.0 ±10.7), are found to impact ILUC-GHG values as these determine how much palm oil is substituted in the international market in response to changes in prices, thereby altering where oilseed production expands at the margin. Land expansion into natural vegetation (35.3 ±10.4) was found to impact ILUC-GHG value ranges by affecting land use expansion dynamics mainly in the US and SEA. The rest of parameters assessed have minor effects in the comparativestatic ILUC-GHG value, mainly determining cropland area requirements to meet oilseed demand for biodiesel uses and further crop demand for animal feed (yield elasticity,demand elasticity for animal products, and exogenous yield projection); as well as the GHG emissions intensity of forestland cover loss and peatland oxidation. The ranges of recursive-dynamic ILUC-GHG values widen for most parameters, particularly, the economic parameters which govern adjustments in international trade and oilseed markets. The trade elasticity of vegetable oils (44.2 ±18.8 gCO2e/MJ) shows the largest variations. The role of parameters determining the magnitude of the livestock rebound effect intensifies with the recursive-dynamic ILUCGHG value. For instance, the yield elasticity (44.7 ±16.1), and the demand elasticity for animal products (44.0 ±14.1) have sizable impacts on the ILUC-GHG value ranges. The demand elasticity for vegetable oils (34.2 ±14.1) and the land expansion into natural vegetation (53.5 ±18.0) remain very influential, while the substitution elasticity among vegetable oils (32.1 ±3.5) matters less for ILUC-GHG value variability relative to the comparative-static ILUC-GHG value. 3.2.2. Monte Carlo Analysis. To test sensitivity of the GLOBIOM results to different assumptions, the OAT estimates presented above were complemented with a full MC analysis targeting the same 11 parameters in the model. Under a MC analysis, each parameter is assigned an assumed Figure 2. Top and bottom panels show the distribution of ILUC-GHG values (gCO2e/MJ) for each parameter in the OAT analysis for the comparative-static and recursive-dynamic ILUC-GHG values. Bars represent the 10th and 90th percentile. Whiskers represent the minimum and maximum values. Solid vertical lines represent the ILUC-GHG value in the central scenario. Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX E probability distribution and a large number of simulations is performed, based on a number of simultaneous randomized draws from each distribution of parameters. This method allows testing of the overall sensitivity of the model results to parametric uncertainty. However, this approach does not touch upon the decision uncertainty related to model design. Therefore, these results do not represent a comprehensive probabilistic estimate of ILUC-GHG values for this modeling framework, but rather a probabilistic estimate specific to the current model structure. While the OAT identifies the effect of varying individual parameters on ILUC-GHG values, the MC analysis allows derivation of uncertainty ranges when all tested Figure 3. Distribution of comparative-static and recursive-dynamic ILUC-GHG values (gCO2e/MJ) and underlying emissions sources obtained from the MC analysis. Boxes show values between the 10th and 90th percentiles. The upper whisker is the maximum and the lower whisker is the minimum value when excluding outliers according to the “1.5 rule;” an estimate is considered an outlier if it is < Q1−1.5 x IQR or > Q3 + 1.5 x IQR, where IQR is the interquartile range. This is meant to exclude results obtained from the combination of the most extreme values of the parameters, given the linear programming formulation of the GLOBIOM model. Figure 4. Distribution of regional comparative-static (upper panels) and recursive-dynamic (lower panels) ILUC-GHG values (gCO2e/MJ) obtained from the MC analysis (left) and decomposition of the corresponding standard deviation by emissions source (right). The boxes in the box and whiskers plot show values between the 10th and 90th percentiles; the upper whisker is the maximum and the lower whisker is the minimum value when excluding outliers according to the “1.5 rule.” An estimate is considered an outlier if it is < Q1−1.5 x IQR or > Q3 + 1.5 x IQR, where IQR is the interquartile range. This is meant to exclude results obtained from the combination of the most extreme values of the parameters, given the linear programming formulation of the GLOBIOM model. South America includes Brazil, Argentina, and rest of South America. Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX F parameters are randomly and simultaneously varied. The MC analysis finds ILUC-GHG values for US soybean biodiesel between the 10th and 90th percentiles ranging from 15.1 to 67.7 gCO2e/MJ in comparative-static and from 8.4 to 77.4 gCO2e/MJ in recursive-dynamic approaches (Figure 3). The full comparative-static (recursive-dynamic) ILUC-GHG value range spans from a minimum value of −17.0 (−23.7) gCO2e/ MJ to a maximum value of 98.7 (112.8) gCO2e/MJ when excluding outliers, with both distributions generally symmetric. The mean values and spread of values slightly increase in the recursive-dynamic ILUC-GHG value (42.4 ±25.9 gCO2e/MJ) relative to the comparative-static one (40.8 ±20.5 gCO2e/ MJ), mainly through increased emissions (and associated uncertainty) from peatland oxidation (Figure 3). Most ILUCGHG values fall on the positive side of the distribution, with 98.8% and 94.7% of runs being above zero in the comparativestatic and recursive-dynamic ILUC-GHG values, respectively. Emissions from natural land conversion and peatland oxidation are the most influential emissions sources, showing a wider distribution than all others (Figure 3), with a standard deviation of 13.5 and 11.4 gCO2e/MJ, respectively, in the comparative-static approach; 18.2 and 12.3 gCO2e/MJ in the recursive-dynamic one. Natural land conversion contributes 47% and 58%; and peatland emissions contribute 34% and 26% of the total variance of the comparative-static and recursive-dynamic ILUC-GHG values, respectively. Despite the large variation, the distributions for these two sources of emissions remain positive within the 10th and 90th percentiles, which is consistent with findings from the central case and OAT analyses (see Sections 3.1 and 3.2.1). Conversely, the distribution of agricultural biomass emissions falls entirely to the left of the y-axis in Figure 3, indicating net global carbon sequestration in agricultural biomass. This stems from a significant portion of the added agricultural biomass coming from palm plantations. Agricultural biomass is rather narrowly distributed around the mean value with a contribution to total variance of 7% (6%) in the comparative-static (recursivedynamic) setup. The remaining emissions from SOC and natural land reversion span from negative to positive and are centered around zero. SOC contributes 12% (10%) of the total variance in the comparative-static (recursive-dynamic) ILUCGHG values, while natural land reversion plays a marginal role (<1% of the variance). As for total ILUC-GHG values, all emissions sources’ distributions are more spread in the recursive-dynamic setup than in the comparative-static one, including SOC (0.3 ±7.6 vs 2.7 ±6.7 gCO2e/MJ), natural land conversion (26.2 ±18.2 vs 27.9 ±13.5 gCO2e/MJ), peatland oxidation (37.4 ±12.3 vs 34.3 ±11.4 gCO2e/MJ), and agricultural biomass (−22.8 ±5.7 vs −21.9 ±5.1 gCO2e/ MJ). This result is due to the increased uncertainty in the two main mechanisms driving ILUC-GHG values, namely the substitution in vegetable oil and feed feedstock markets, as well as in the subsequent livestock rebound effect. The widest uncertainty ranges in both comparative-static and recursive-dynamic ILUC-GHG values, respectively, are found for SEA (38.0 ±15.3; 41.8 ±18.1 gCO2e/MJ) and SAM (−7.2 ±7.4; −5.3 ±14.1 gCO2e/MJ), while mean ILUC-GHG values for the US remain rather low and narrow (4.7 ±3.7; 6.4 ±4.4 gCO2e/MJ) (Figure 4). SEA respectively contributes 73.5% and 54.8% of the comparative-static and recursive-dynamic ILUC-GHG value variance. Emissions from peatland oxidation are again a key driver of the variation, contributing 54% (48%) of the variance of SEA’s comparativestatic (recursive-dynamic) ILUC-GHG value, followed by natural land conversion with 25% (34%). For SAM, regional ILUC-GHG value variation is less pronounced than for SEA (Figure 4). ILUC-GHG values for SAM are negative in 71% (60%) of the MC simulation runs in the comparative-static (recursive-dynamic) setup. These negative ILUC-GHG values correspond to those MC simulations with lower deforestation rates in Brazil and Argentina, coupled with higher assumed EF from forest biomass loss values, which translates into net carbon gains compared to the baseline. MC runs that result in cropland expansions in SAM are largely showing higher demand elasticities for vegetable oils, higher yield projection responses, lower vegetable oil substitution elasticities, and lower yield elasticities (Section S2.4). This leads to a marginally greater demand for soybean, corn, and soybean Figure 5. Marginal effects on the ILUC-GHG value (gCO2e/MJ) from the variation of key model parameters in the MC analysis. Only the most influential parameters are shown. The bands around the lines indicate the 95% confidence interval of the marginal effect. The values in the x-axis indicate the percentile within the defined distributions of shifter values applied to the parameters as described in Table S1. Agricultural biomass emissions factor: EF of agricultural biomass in palm plantations; forest emissions factor: EF from forest biomass loss; palm to peat expansion response: expansion response of palm into peatland. Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX G oil from SAM used to meet demand for food, feed, and other uses in regions outside of the US, while palm oil continues to have a limited market penetration in the global mix. The BART analysis shows that the EF from forest biomass loss causes the widest marginal variation of ILUC-GHG values under the comparative-static approach with a 95% confidence interval; while that very parameter and especially the expansion response of palm into peatland cause the sharpest variation in recursive-dynamic ILUC-GHG values, as the latter directly influences uncertainty in peatland emissions. The marginal effects of the EF of agricultural biomass in palm plantations and the peatland EF also become marked under the two approaches (Figure 5). The yield projection response has a more irregular effect, especially under the comparative-static approach, given the role that this parameter plays in determining cropland area dynamics in US and SAM. On the one hand, the highest shifter values result in net cropland increases in SAM; on the other hand, the lowest values imply larger areas converted into soy in US. Under the recursive-dynamic approach, yields tend to increase for all regions and crops over the period considered, contributing to mitigate global cropland area expansion and thus ILUC as the shifter increases. The marginal effect of the economic parameters increases under the recursive-dynamic approach, though remains more moderate than that of the biophysical ones. The marginal effect of the vegetable oil substitution elasticity becomes especially uncertain with the highest values. Unlike the OAT analysis, BART underlines the more prominent role of biophysical parameters in driving ILUC results uncertainty when propagated with MC simulations. 4. DISCUSSION 4.1. Major Market-Mediated Responses. Increasing demand for biofuels produced from oilseeds generates diverse market-mediated effects and cross-sectoral interactions, which can ultimately translate into intensified competition for land among all uses and increased land prices. 72,79−81 These characteristics make the potential spillover effects of vegetable oil demand shocks complex, establishing the need for systematic uncertainty analysis when estimating ILUC-GHG values with global economic models. Our study improves the understanding of the responses across vegetable oil and oilseed markets and provides additional evidence that modeled soybean biodiesel ILUC-GHG values are highly sensitive to both economic and biophysical parameters. OAT analysis and the combination of BART and MC analyses help identify model parameters related to vegetable oils markets, and in particular GHG emissions from oil palm expansion, as key determinants of ILUC-GHG values and ranges. Results highlight the role of particular regions, SEA and SAM, where land use spillovers from the US soybean oil biodiesel consumption shock are significant. The level of estimated cropland expansion and associated pressures on natural ecosystems in SEA and SAM in modeled results depends significantly on the assumed vegetable oil substitution elasticity − the ease with which markets may shift consumption from one vegetable oil to another in response to demand shocks−and the vegetable oil demand elasticity −the ease with which consumption of vegetable oils overall may decrease when prices rise. These parameters alter how much additional soybean and palm oil is produced in SAM and SEA respectively, thereby contributing to the interplay between vegetable oil demand in the US and enhanced production abroad. Across the MC runs, GLOBIOM consistently estimates increasing palm oil use as the primary substitute for diverted soybean oil from food and other uses in international markets, causing natural land conversion and peatland emissions in SEA. This is consistent with several other CLCA studies that identify palm oil from SEA as the most cost-competitive but emission intensive 82 marginal source of vegetable oil feedstock in the market. 19,45,83,84 This market-mediated effect is also identified in GTAP-BIO studies, 55 especially in scenarios with relatively higher substitution elasticities that only consider soy oil-palm oil substitution. The share of Malaysia and Indonesia in the estimated ILUC-GHG emissions value for soybean biodiesel (17.5 g CO2e/MJ, + 2 billion gallons) is 78%; in our case it is <50% (see Figure 1). If other vegetable oils and animal fats are taken into account, most of the additional demand for soybean oil is diverted to those produced in the US or in regions other than the SEA. 55 This highlights the role of the elasticity of substitution among the various types of biodiesel feedstocks, especially when other kinds of vegetable oils and animal fats are considered. In this way, we could expect lower ILUC-GHG values when including more feedstocks available in the US, such as tallow or used cooking oil, as combined with greater assumed elasticities. This would allow for nonpalm vegetable oils and animal fats to substitute for soybean oil more than in our results, hence decreasing demand for palm oil and land conversion in SEA. This difference among findings highlights the decision uncertainty existing among models, which may be partially related to the elasticity of substitution among the various types of vegetable oils and animal fats in the US; greater assumed elasticity may allow for nonpalm vegetable oils and animal fats to substitute for soybean oil more so than we have found in these results. Conversely, GLOBIOM generally estimates a slowing rate of natural land conversion in SAM. Brazil loses market share in global soybean markets through the reduced demand for Brazilian soybean meal for feed applications due to the increased availability of relatively cheap US soybean meal. The MC analysis still finds a minority of runs (6.3% in 2030, and 45.9% in 2050), in which cropland expands in SAM, increasing ILUC-GHG values. This effect emerged across diverse arrays of input parameter values and could not be explained by any one specific parameter. Since both SAM and SEA include carbon-rich forests and natural lands, any cropland expansion in these regions will impact the ILUCGHG value of US soybean biodiesel (Figure S7) and associated uncertainty ranges (Figure 4). This highlights the importance of the four biophysical parameters for the MC analysis, as these determine the net GHG emissions per hectare when forest is converted into cropland and the total peatland emissions from palm expansion. GLOBIOM consistently estimates additional cropland expansion in major grain-producing regions outside the US to meet increased demand for animal products when soybean meal availability increases. This livestock rebound effect 78 is also subject to uncertainty in the parameters considered and contributes additional GHG emissions in regions such as SAM, SAS, and the US, especially under the recursive-dynamic approach in which the impacts of uncertainty in these parameters are accounted for over multiple decades. Methodological Considerations. The market-mediated responses summarized above highlight the integration of Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX H globalized food, feed, and fuel markets. Similar effects could be expected when simulating increased demand for other oilseeds, especially those with meal coproduction, such as rapeseed. 40 Uncertainties persist in estimates of projected future impacts from vegetable oil-based fuels. Our results highlight the need for better empirical estimates for critical economic parameters such as elasticities of vegetable oil substitution and consumer demand elasticities to vegetable oil prices. Further efforts are needed to refine estimates on biophysical parameters, such as EFs from land conversion benefiting from improving spatially explicit data on yields, land areas, and associated carbon stocks. 31 Increased availability of remote sensing data and associated processing capacities can help to better characterize the inherent uncertainty related to the variability of the modeled systems, e.g., spatial variability in carbon stocks. 25 , 42 However, given the size and mathematical complexity of global economic models, epistemic uncertainty will always remain a challenge when estimating and interpreting ILUC-GHG values. MC analysis is often used to assess sensitivity and uncertainty in results from global economic models. However, the number of varied model parameters tends to be limited. This study uses MC analysis to explore the ILUC-GHG value sensitivity to parametric uncertainty, where both the number of parameters and their individual value ranges and distributions have been chosen in advance, essentially defining the assumed likelihood that each parameter takes a value within a certain range. Although this is a common practice for uncertainty analysis, 14,49,85 other important parameters may be overlooked by this approach. It should also be noted that for many parameters, the probability distribution is not very well-known, such that epistemic uncertainty may not be fully represented. The decision of which parameters to include and how to shape their stochastic distribution of values is a critical analytical choice by the modeler. This highlights the importance of transparency regarding the choice of parameters, model structure and underlying assumptions for CLCA, while assessing trade-offs on different temporal and spatial scales. 37,42,43,86 The MC analysis presented helps understand model behavior when specific model inputs are modified arbitrarily, as a mean to propagate uncertainty in ILUC-GHG values. However, a more comprehensive uncertainty analysis would require a previous characterization of an expected range of possible behaviors or data and would need to carry that characterization through to a range of implied possible outcomes. 64,87 Despite the considerations above, this study covers a comprehensive set of parameters and distributions in the MC simulation based on recent literature and expert judgment (Table S1 in ESM). Seven economic parameters and four biophysical parameters are varied, expanding the number of parameters varied relative to previous GLOBIOM uncertainty studies. 40 The assumed distributions could be refined further, for instance, by defining the elasticities based on econometric methods and empiric approaches, or by using improved, finerscale, and more up-to-date remote sensing data to determine the stochastic distribution of biophysical parameters across the globe. This is out of the scope of the present study but should be addressed in future work. Given the nonlinear nature of ILUC-GHG values, our MC analysis shows that including more parameters does not necessarily yield wider ILUC-GHG value ranges, as many parameters interact, especially the economic ones that affect several regions and commodities simultaneously. This study applies BART as a tool to better understand the role of each parameter in the variation of ILUC-GHG values, aiding the interpretation of MC analysis results from global economic models. Still, MC analysis proves useful to propagate uncertainty in input parameters and obtain probability distributions of ILUC and other model-based outcomes. 86,88 Other crucial modeling choices include the magnitude of feedstock consumption simulated exogenously, the biofuel shock size, and the baseline used as counterfactual. 18,40 Furthermore, the time frame over which a biofuel shock is evaluated can have a significant impact on resulting ILUCGHG value. 86 These choices explain differences in the mean ILUC-GHG value estimated in the MC analysis in this study (between 40.8 and 42.7 gCO2e/MJ, depending on the approach) as compared to previous assessments such as Prussi et al. 48 The comparative-static approach may be more appropriate when one seeks to understand near-term market dynamics, as it relies on only assumptions relevant to the next immediate future time step and provides better computational efficiency. The recursive-dynamic approach may be more appropriate when seeking to understand longer-term bioeconomic dynamics, such as impacts on ILUC-GHG emissions, as it better captures the inherent temporal variation and uncertainty in such estimates. Most recent literature has considered only one of these approaches, with the predominant focus on the comparative-static approach. The two metrics apply different temporal scopes, hence providing different insights since results vary depending on the baseline developments and underlying assumptions. For example, the recursive-dynamic ILUC-GHG value (up to 2050) also considers other important drivers over time, e.g., the fact that yields are larger in 2050, which reduces the impact of the biofuel shock. This study’s explicit comparison of comparativestatic and recursive-dynamic approaches in an otherwise consistent modeling framework demonstrates the methodological choices modelers should consider in the face of model uncertainty, particularly related to the treatment of time. When relying on biofuel CLCA modeling to inform decision-making, modelers should acknowledge which uncertainties are and are not represented in their supporting analyses, including related to assumptions implicit to utilized methods of quantifying ILUC-GHG impacts. ■APPENDIX CONV: natural land cover conversion including aboveand below-ground biomass, dead wood, litter, and harvested wood products REV: natural land cover reversion or forest regrowth BIOM: carbon sequestration in agricultural biomass of bioenergy crops SOC: soil organic carbon PEAT: GHG emissions from peatland oxidation = [ + + ] × LUC MtCO CONV MtC BIOM MtC REV MtC ( ) ( ) ( ) ( ) 44 12 y y y y 2 (1) Environmental Science & Technology pubs.acs.org/est Article https://doi.org/10.1021/acs.est.3c09944 Environ. Sci. Technol. XXXX, XXX, XXX−XXX I