Global cooling induced by biophysical effects of bioenergy crop cultivation

Bioenergy crop with carbon capture and storage (BECCS) is a key negative emission technology to meet carbon neutrality. However, the biophysical effects of widespread bioenergy crop cultivation on temperature remain unclear. Here, using a coupled atmosphere-land model with an explicit representation of lignocellulosic bioenergy crops, we find that after 50 years of large-scale bioenergy crop cultivation following plausible scenarios, global air temperature decreases by 0.03~0.08 °C, with strong regional contrasts and interannual variability. Over the cultivated regions, woody crops induce stronger cooling effects than herbaceous crops due to larger evapotranspiration rates and smaller aerodynamic resistance. At the continental scale, air temperature changes are not linearly proportional to the cultivation area. Sensitivity tests show that the temperature change is robust for eucalypt but more uncertain for switchgrass among different cultivation maps. Our study calls for new metrics to take the biophysical effects into account when assessing the climate mitigation capacity of BECCS.

L arge-scale bioenergy crop cultivation with carbon capture and storage (BECCS) has been identified by integrated assessment models (IAMs) as a major negative emission technology (NET) for removing CO 2 from the atmosphere 1-3 . Most mitigation scenarios which limit the global temperature increase below 2 or 1.5°C in 2100 assume widespread deployment of this NET 2,4 . In the scenarios where BECCS was the only land-based negative emission option 5,6 , about 128 PgC needs to be removed through BECCS to reach the temperature limiting target under Representative Concentration Pathway 2.6 (RCP 2.6). However, the real mitigation potential of BECCS depends on multiple dimensions, such as achievable bioenergy crop yields 7 , available cultivation lands 8 , land use change emissions 5,6 , and supply of water 9,10 and nutrients [11][12][13] . Marginal lands (e.g., abandoned agricultural land or degraded lands) are targeted for bioenergy crop cultivation to avoid land competition with food crops and forests 8,14,15 , but large-scale use of current agricultural systems for bioenergy crops was also proposed 12,16,17 . In addition, the biophysical impacts of bioenergy crops are important to the mitigation potential of BECCS. Biophysically cooling effects were found in the US when converting annual crops to perennial bioenergy crops 18 , but the biophysical interactions between largescale bioenergy crop cultivation and climate remain unclear at the global scale. Biophysical climate effects could either partly offset or enhance the climate benefits from carbon-dioxide removal (CDR). Changes in land management and land cover alter land surface properties like albedo and aerodynamic resistance and influence the surface energy budget [19][20][21][22][23] . In turn, these changes modify air temperature (T a ) locally and regionally, in the latter case through atmospheric circulation changes and teleconnection mechanisms 24 . It is thus important to quantify the biophysical climate effects of large-scale bioenergy crops to fully assess their role in climate mitigation.
Here, we simulated the biophysical climate impact of a range of future bioenergy crop cultivation scenarios (Supplementary Methods 1) using an Earth system model (IPSL-CM 25 ). IPSL-CM deploys a land surface scheme with an explicit representation of the main lignocellulosic bioenergy crops (ORCHIDEE-MICT-BIOENERGY 26 ) as four plant functional types (PFTs): eucalypt, poplar & willow (sharing the same PFT because of their similar ecological processes and parameters), miscanthus, and switchgrass 26 . The resolved management practices include the periodic harvest of woody crops, separate harvested biomass pools, and regrowth from belowground biomass after harvest. The parameters controlling physiological processes, growth, and phenology were calibrated using plant measurements and a global database of yields from field trials 26,27 . Model performance was evaluated against those field measured yields 27 and observationbased biophysical variables, such as albedo and evapotranspiration, which control the energy budget (see Supplementary Methods 2, Supplementary Tables 1 and 2, and Supplementary Figs. [1][2][3][4][5][6][7][8][9]. From the coupled simulations, we find that large-scale bioenergy crop cultivation induces a biophysical cooling effect at the global scale, but the air temperature change (ΔT a ) has strong spatial variations and interannual variability. Compared to the herbaceous crops, changes in the energy fluxes induced by woody crops in the cultivation regions are larger, and the cooling effect is more robust across different cultivation maps.

Results
Spatial patterns of biophysical ΔT a . A composite map (Supplementary Fig. 10f) for bioenergy crop deployment around the world was generated from three marginal land datasets 8 Supplementary Fig. 11), and four idealized bioenergy crop scenarios (i.e., S euc , S p&w , S mis , and S swi ) where all BECCS grid cells are assumed to be fully covered by a single type of the four bioenergy crops (i.e., eucalypt, poplar & willow, miscanthus, and switchgrass), with the remaining land cover being the same as in S ref . Differences (Δ) in simulated climate between each bioenergy scenario and the reference scenario represent the biophysical climate effects of each bioenergy crop (see "Methods" and Supplementary Methods 4.1).
Although the assumed bioenergy crop cultivation accounts for 3.8% of global land area in the composite cultivation map (Fig. 1a), air temperature is significantly (p < 0.1) changed over about 10% of the global land area (grid cells with black shading in Fig. 1c-f, 84-89% of which are distributed outside of the BECCS regions). At the global scale, different bioenergy crop cultivations result in cooling effects with ΔT a = −0.03, −0.08, −0.06, and −0.04°C for S euc , S p&w , S mis , and S swi , respectively (Figs. 1c-f and 2a), when averaged over the final 10 years of simulations. The relatively small magnitude of global mean ΔT a results from the contrasting warming and cooling effects in different grid cells ( Fig. 1c-f). In fact, the spatial variations (standard deviation, SD) of ΔT a range from 0.38 to 0.42°C across different scenarios, much higher than the global mean ΔT a . Strong interannual variability (IAV) was also found in the last 20 years of the global mean ΔT a (IAV = 0. 16 Spatial patterns of ΔT a differ among the four idealized bioenergy crop scenarios using the composite cultivation map (Fig. 1c-f). T a changes more in cold regions than in warm regions ( Fig. 1b and Supplementary Fig. 21). In boreal regions, all bioenergy crops exert a cooling effect on average (from 50°N to 80°N in Fig. 1b), mainly due to strong cooling signals in northern Europe and Siberia ( Fig. 1c-f). In the warm northern temperate regions, the impacts of bioenergy crop plantation scenarios on T a are very different (Fig. 1b). We find a cooling effect in the band from 30°N to 50°N for S p&w , with cooling epicenters in China, central Asia and the eastern US. In contrast, when switchgrass or eucalypt plantations are assumed (i.e., S swi and S euc ), we obtained warming effects in these regions (Fig. 1c, f). Across temperate Asia, the warming impact in S euc (Fig. 1c) is stronger than in scenarios with herbaceous bioenergy crops (i.e., S mis and S swi , Fig. 1e, f). In tropical regions, the magnitude of ΔT a is relatively small for the scenarios based on the composite cultivation map (Fig. 1b). Compared to the herbaceous bioenergy crop scenarios (i.e., S mis and S swi ), the woody bioenergy crop scenarios (i.e., S euc and S p&w ) significantly cool the climate in regions of widespread cultivation, e.g., in southeast Brazil (Fig. 1c, d). Similarly, in sub-Saharan Africa, woody bioenergy crop scenarios have significant cooling signals, while temperature change from herbaceous bioenergy crops is weaker (Fig. 1c-f). Bioenergy crops have little impact on T a in southeast Asia. A cooling effect is found in the north of Australia for all bioenergy scenarios using the composite cultivation map, while central and west Australia and New Zealand generally show warming effects, except in S mis (Fig. 1c, d, f). By comparison, S mis induces more cooling effects than warming effects in Oceania (Fig. 1e).
Contribution of energy fluxes to ΔT a . We decomposed the ΔT a induced by bioenergy crops into components related to the local surface energy balance and to the large-scale circulation 23,24 , based on: where ΔT local a is the temperature change induced by changes in the local surface energy balance (see details in Eqs. (13) and (14) of "Methods" and Supplementary Figs. 14-18). ΔT cir a represents the impact of atmospheric circulation changes induced by changing the coverage of bioenergy crops (see details in "Methods" and Supplementary Methods 2.4). ΔT local a can be determined from refs. 23,24 : where f a is an energy redistribution factor (see "Methods"), and ΔSW dn , Δα, ΔLE, ΔG, Δr a , and Δε a are changes in surface downward shortwave radiation, surface albedo, latent heat flux, ground heat flux, aerodynamic resistance, and longwave radiative emissivity of air, respectively. C p , ρ, ε s , and σ are specific heat capacity of air, air density, surface longwave radiative emissivity (prescribed as 1.0 in the coupled model), and Stephan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ), respectively.
Globally, for woody bioenergy crop scenarios based on the composite cultivation map (i.e., S euc and S p&w ), ΔT local a is a net cooling effect with dominant contributions from decreased aerodynamic resistance (S euc ) or increased albedo (S p&w ). On the other hand, the sign of ΔT cir a differs between S euc (warming) and S p&w (cooling) (Fig. 2a). Compared to S p&w , the cooling effect of ΔT local a in S euc is largely counteracted by the warming effect from ΔT cir a , leading to a smaller net ΔT a . In contrast, for herbaceous bioenergy crop scenarios, ΔT local a is a net warming effect, while ΔT cir a is a cooling effect (Fig. 2a). For both S mis and S swi , the cooling effect of ΔT cir a outweighs the warming effect of ΔT local a causing a global net cooling effect (global mean ΔT a < 0). The warming effect of ΔT local a for the herbaceous bioenergy crop scenarios is mainly contributed by increased aerodynamic resistance, decreased latent heat and increased downward shortwave radiation. The global average cooling effect of albedo and ground heat flux and the warming effect of downward shortwave radiation are consistent among the four bioenergy crop scenarios based on the composite cultivation map (Fig. 2a).
Over the cultivated regions in the composite cultivation map, bioenergy crop cultivations show a cooling effect except for switchgrass (Fig. 2b). S euc , S p&w and S mis enhance local evapotranspiration and thus decrease T local a ( Fig. 2b and Supplementary Figs. 14-17), whereas S swi decreases evapotranspiration which offsets the cooling effects from other terms in Eq. (2) (e.g., a higher albedo), resulting in an increased ΔT local a ( Fig. 2b and Supplementary Fig. 15). All the idealized bioenergy crop scenarios using the composite cultivation map reduce the aerodynamic resistance which induces a cooling effect where they are cultivated. In addition, the magnitudes of the different gross components in Eq. (2) are much larger for woody bioenergy crops than for herbaceous crops over the cultivated regions ( Fig. 2b), due to the differences in vegetation properties (e.g., LAI and roughness) and changes in the related variables (e.g., wind speed, humidity, and cloud) between woody and herbaceous crops (Supplementary Figs. 22 and 23, see details in Supplementary Discussion 1).
Outside the regions of bioenergy crop cultivation, ΔT local a shows warming effects for all the idealized bioenergy crop scenarios based on the composite cultivation map (Fig. 2c), while ΔT cir a has cooling effects. Atmospheric circulation redistributes surface energy changes 23,28,29 , and wind speed and surface roughness height are tightly associated with atmospheric circulation ( Supplementary Fig. 24). Therefore, outside the cultivated regions in the composite cultivation map, where there is no land cover change in the bioenergy crop scenarios (same as S ref ), atmospheric circulation and teleconnections alter local radiative forcing indirectly through exchanges of heat between different regions (Supplementary Figs. [14][15][16][17]. Changes in local radiative forcing further influence vegetation gas exchange and structure and thus the surface energy balance, e.g., changes in downward shortwave radiation, air emissivity and aerodynamic resistance ( Fig. 2c and Supplementary Figs. 12-15). This is supported by the fact that LAI and radiative forcing are slightly changed outside the regions of bioenergy crop cultivation (Supplementary Fig. 24 and Fig. 2c).
T a changes at continental scale. Based on the composite bioenergy crop cultivation map, the effects of local surface energy balance changes in air temperature (ΔT local a ) range from warming in Oceania, which has the smallest bioenergy crop cultivation area, to  , and outside the BECCS regions (c). In each panel, four groups of bars represent four idealized bioenergy crop scenarios (i.e., eucalypt, poplar & willow, miscanthus, and switchgrass) based on the composite cultivation map. Bars for ΔT cir a , G, ε a , r a , LE, SW dn , and α represent different components that contribute to the air temperature change, i.e., the altered atmospheric circulation, ground heat flux, longwave radiative emissivity of air, aerodynamic resistance, latent heat flux, surface downward shortwave radiation, and surface albedo, respectively. The net air temperature change is shown by red dots, which is the sum of all bars. Black asterisks represent the contribution of altered local surface energy balance, which is the sum of the colored bars excluding the gray bar. Note that the scale of the x-axis is different for the three panels.
cooling in Europe, which has the largest area of bioenergy crop cultivation (Fig. 3a). In fact, ΔT local a is significantly correlated with the area of bioenergy crop cultivation in the composite cultivation map across different continents (p < 0.05) for woody bioenergy crop scenarios ( Supplementary Fig. 25). However, changes in atmospheric circulation (ΔT cir a , Fig. 3b) that redistribute energy spatially show a more subtle pattern (Fig. 3b), either enhancing or weakening the ΔT local a signal. As a consequence, mean ΔT a at the continental scale depends on both ΔT local a and ΔT cir a as well as on the bioenergy crop type considered. Therefore, continents with a larger area of bioenergy crop cultivation do not always show the stronger cooling effect ( Fig. 3c and Supplementary Fig. 25).
Woody bioenergy crop cultivation largely decreases T local a in Europe (ΔT local a = −0.16 and −0.12°C for S euc and S p&w , Fig. 3a and Supplementary Fig. 26) due to the effects of decreased aerodynamic resistance, increased evapotranspiration, and decreased downward shortwave radiation (Supplementary Fig. 26). The continental mean ΔT a is also contributed by ΔT cir a (0.06 and −0.06°C for S euc and S p&w , Fig. 3b and Supplementary Fig. 26). In South America, which has substantial areas of bioenergy crop cultivation in our scenarios (133 M ha), the cooling effect of ΔT local a is counteracted by a warming effect of ΔT cir a , leading to a net increase of T a in all scenarios based on the composite cultivation map except S euc (Fig. 3 and Supplementary  Fig. 26). Bioenergy crop scenarios generally result in weak effects on the mean T a of Africa, with all components contributing little to ΔT a (Supplementary Fig. 26). Asia and North America have similar bioenergy crop cultivation areas (40 and 41 M ha, respectively) and both show net cooling effects for all bioenergy crop scenarios, but the signs of ΔT local a and ΔT cir a are generally opposite in these two continents ( Fig. 3 and Supplementary Fig. 26). In contrast, in Oceania, which has the smallest bioenergy crop cultivation areas, all bioenergy crop scenarios using the composite cultivation map, with the exception of S mis , show warming effects (ΔT a > 0) contributed by positive ΔT cir a (Supplementary Fig. 26).  In addition, although mean T a over the regions of bioenergy crop cultivation based on the composite cultivation map within each continent generally decreases due to reduced aerodynamic resistance and enhanced latent heat, especially for woody bioenergy crops ( Supplementary Fig. 27), the continental mean ΔT a is mainly contributed by ΔT a outside the cultivated regions, implying the importance of the teleconnection of temperature changes ( Supplementary Fig. 28). For example, poplar & willow cultivation decreases the mean T a over their cultivated regions in all continents (Supplementary Fig. 27), but this signal is offset by teleconnection warming outside the cultivated regions in Oceania, Africa, and South America ( Supplementary Fig. 28), resulting in a positive mean ΔT a for the whole continent ( Fig. 3c and Supplementary Fig. 26). Despite contrasting changes in T a within and outside the bioenergy crop cultivation regions in the composite cultivation map, the mean ΔT a and the contributions of various energy fluxes for the whole continent ( Supplementary  Fig. 26) generally agree with the changes in the grid cells with significant ΔT a (grid cells with black shading in Fig. 1c-f) Supplementary Fig. 10a, c, d) and two representative bioenergy crop types (eucalypt and switchgrass). In the additional simulations based on the three individual cultivation maps and two bioenergy crop types, the global net ΔT a (averaged over the last 10 years) ranges from −0.05 to 0.05°C across various simulations with significant ΔT a appearing in 9-25% of the global lands ( Fig. 4 and Supplementary Fig. 30). Similar to the results from the composite map (Fig. 1), ΔT a also has strong spatial heterogeneity (spatial SD ranges from 0.25 to 0.37, Supplementary Fig. 30) and high interannual variability (20-year IAV ranges from 0.15 to 0.2°C across different scenarios).
At the global scale, the signs of ΔT local a and most energy components for eucalypt are robust among all cultivation maps, but the magnitudes vary, especially for ΔT cir a (Fig. 4). Different from the net cooling effect of ΔT a using composite, IMAGE, or MAgPIE map, the positive ΔT a in S Cam euc (i.e., the simulation using the Campbell-high map for eucalypt) is caused by the large warming effect of ΔT cir a , which offsets the cooling effect of ΔT local a . By contrast, the changes of biophysical effects of switchgrass are different among simulations with different cultivation maps. ΔT a shows a cooling effect in S swi , but a warming effect in S IMA swi and S MAg swi , and a nearly neutral effect in S Cam swi . Compared to eucalypt, the responses of each component to switchgrass cultivations are relatively small using all cultivation maps, resulting in higher sensitivity of ΔT a to cultivation maps. The signs of ΔT a for switchgrass are largely determined by ΔT cir a . Similarly, at the continental scale, the signal of ΔT local a for eucalypt is generally robust with respect to the choice of a cultivation map, showing consistent patterns between the idealized composite maps and the three individual maps in most continents, but the signs of ΔT a also depend on the ΔT cir a (Supplementary Fig. 31). In agreement with the differences between eucalypt and switchgrass at the global scale (Fig. 4), the magnitudes of biophysical effect change for switchgrass are also lower than eucalypt at the continental scale, and the patterns are more diverse across various maps (Supplementary Fig. 31 and 32). Only ΔT local a in Europe and ΔT cir a in Oceania and Africa are consistent for switchgrass cultivation across the four maps (one composite and three individual maps). Combining ΔT local a and ΔT cir a , however, no consistent ΔT a for switchgrass were found in any continent ( Supplementary Fig. 32).
The original land use types upon which bioenergy crop cultivation is expanded (Supplementary Fig. 10) also determine the changes in the surface energy budget. Cultivation regions from the MAgPIE and Campbell-high maps were generally converted from short vegetation like cropland and grassland, while in the IMAGE map, 78% of the bioenergy cultivation lands were converted from forest. Even so, both ΔT local a and ΔT a in S IMA euc are negative (Fig. 4), indicating the stronger cooling effects of eucalypt plantations than previous forest types. Compared to S euc Fig. 4 Contributions of different components to air temperature changes at the global scale using different cultivation maps. Two panels represent eucalypt cultivations (a) and switchgrass cultivation (b). Four rows of bars represent four cultivation maps (i.e., the composite, MAgPIE, IMAGE, and Campbell-high map). Bars forΔT cir a , G, ε a , r a , LE, SW dn , and α represent different components that contribute to the air temperature change, i.e., the altered atmospheric circulation, ground heat flux, longwave radiative emissivity of air, aerodynamic resistance, latent heat flux, surface downward shortwave radiation, and surface albedo, respectively. The net air temperature change is shown by red dots, which is the sum of all bars. Black asterisks represent the contribution of altered local surface energy balance, which is the sum of the colored bars excluding the gray bar. and S MAg euc , where eucalypt was cultivated on short-vegetation lands, the cooling effects in S IMA euc is even stronger (Fig. 4), probably due to the larger cultivation area in the IMAGE map (523 M ha) than the composite and MAgPIE maps (466 and 408 M ha, Supplementary Fig. 10). By contrast, replacing forests with switchgrass in S IMA swi , both ΔT local a and ΔT a show warming effects (ΔT local a = 0.01°C, ΔT a = 0.05°C). The Campbell-high map represents a widely spreading BECCS cultivation scenario, and its biophysical effects are rather different from those with cultivation areas more concentrated in a few regions (e.g., the larger ΔT cir a in the eucalypt cultivation scenario based on the Campbell-high map, see Fig. 4a and Supplementary Discussion 2). This illustrates the importance of the spatial cultivation patterns on the biophysical effects. Note the differences in biophysical effects between the composite map and the three individual maps also include the differences between the two reference simulations (see "Methods" and Supplementary Discussion 3). For example, the large difference of ΔT a between the three individual maps and the idealized composite map in the boreal regions (Supplementary Fig. 30) is mainly caused by the difference between the two reference simulations (Supplementary Fig. 33).

Discussion
We found that, even though the cultivation area occupies 3.8% ± 0.5% of the global total land area, bioenergy crops exert strong regional biophysical effects, leading to a global net change in air temperature of −0.08~+0.05°C. The CDR target of the BECCS scenarios used in IAMs is 128 PgC till 2100 for RCP 2.6 6 , approximately corresponding to a cooling effect of 0.09~0.32°C (Supplementary Table 5). Thus, the biophysical cooling or warming effects of bioenergy crop cultivation can significantly strengthen or weaken the effectiveness of BECCS in limiting the temperature increments depending on the cultivation map and the bioenergy crop type. However, the temperature changes in the bioenergy crop scenarios do have very large spatial variations and important climate teleconnections to other areas of the globe (Fig. 1c-f). For example, from the four idealized bioenergy crop scenarios based on the composited cultivation map, warming effects in Alaska and northwestern Canada may cause greenhouse gas release from thawing permafrost 30,31 , while strong cooling effects in Eurasia, between 60°N and 80°N, may protect permafrost from thawing or reduce methane emissions from wetlands. Therefore, the biophysical effects of bioenergy crop cultivation do not only alter global temperature directly, but also induce secondary effects on natural greenhouse gas fluxes, which should also be taken into account when considering large-scale BECCS deployment. Meanwhile, altered atmospheric circulation is a major contributor to the net change in global air temperature, but its spatial variations may rely on the choice of the atmosphere model, and the relevant uncertainties need further investigation.
In addition, the six additional scenarios using three plausible cultivation maps and two representative bioenergy crops demonstrated the importance of the crop type choice, the original land use type upon which bioenergy crop are expanded, the total cultivation area and its spatial distribution patterns. Cultivating eucalypt shows generally cooling effects from ΔT a and ΔT local a (Fig. 4), which are more robust than if switchgrass is used as the main bioenergy crop across various cultivation maps, implying that eucalypt is superior to switchgrass in cooling the lands biophysically. Replacing forests with switchgrass in S IMA swi not only results in biophysical warming effects (Fig. 4) but could also release more carbon through deforestation than converting other short vegetation to bioenergy crops. Thus, deforestation should be avoided. The magnitude of changes in the biophysical effects also depends on the total cultivation area. With the largest cultivation areas in the IMAGE map, ΔT a shows greatest cooling effects for eucalypt and greatest warming effects for switchgrass among the four maps (Fig. 4). Last, cultivation concentrated in a few regions rather than spreading fragmented cultivation (like the spatial cultivation patterns in the Campbell-high map) is recommended, which is also conducive to land management.
Air temperature changes in different continents are disproportionate to their areas of bioenergy crop cultivation. In the composite cultivation map, the bioenergy crop cultivation area in Europe (157 M ha) is about four times larger than that in Asia (40 M ha), but both continents have similar cooling magnitudes for all four idealized bioenergy crop scenarios. In other continents, air temperature changes depend on bioenergy crop types. Therefore, there is a requirement for new trans-regional metrics to appreciate the climate impacts of BECCS and for international collaboration to share the responsibility of bioenergy crop cultivation and to address the challenges inherent in slowing down global warming by considering both the biogeochemical and biophysical effects of climate mitigation options. In addition to the biophysical effects, the feasibility of BECCS also relies on other dimensions like the related land use change carbon emissions 6 , the water demand 9,10 and the nutrient cycle [11][12][13] . A full assessment of the nexus of these aspects is thus needed. ORCHIDEE-MICT-BIOENERGY includes specific biomass harvest modules for bioenergy crops: annual harvest of herbaceous crops and the periodic harvest of woody crops after each 5-year rotation. Parameters related to critical processes of bioenergy crops, i.e., photosynthesis, carbon allocation, phenology, and biomass harvest were systematically calibrated based on field data 26,27 . The model is capable of capturing the growth dynamics of bioenergy crops and reproducing the biomass yields in a large global yield observation dataset 26,27 . The LMDz used here is version 6 of the Laboratoire de Météorologie Dynamique atmospheric general circulation model with zooming capability 32,33 . It simulates the fundamental dynamical and physical processes of the atmosphere (e.g., radiation, moist convection, cloud macro-, and micro-physical properties, etc.) using 79 layers in the vertical.
Based on the composited cultivation map, we performed five 50-year coupled simulations starting from 2015: a reference scenario (S ref ) and four idealized bioenergy crop scenarios with cultivation of different bioenergy crops (S euc : eucalypt, S p&w : poplar & willow, S mis : miscanthus, S swi : switchgrass, see Supplementary Methods 4.1). In S ref , the BECCS regions (Fig. 1a, Supplementary Methods 3) were covered by generic grain crops because the BECCS regions in the composite cultivation map were mainly abandoned cropland. In each of the four bioenergy crop scenarios, the BECCS regions in the composite cultivation map were covered only with the corresponding bioenergy crop type. At the end of the 50-year simulations, the key vegetation features (e.g., photosynthesis and leaf area index) that are closely related to the energy budget (e.g., evapotranspiration and albedo) reach a steady state, implying that the simulation length is sufficient to derive a robust signal of biophysical effects from bioenergy crop (see details in Supplementary Methods 4.3 and Supplementary Figs. 12 and 13).
We also performed seven additional simulations to analyze the sensitivity of biophysical effects to different cultivation maps. Seven additional simulations include one reference scenario (S  Table 4). In these three individual cultivation maps, the bioenergy crop coverage is fractional in each grid cell ( Supplementary Fig. 10a, c, d), which is considered more realistic than in the composite map where entire grid cells in the BECCS regions were fully covered by bioenergy crops (Supplementary Fig. 10f). A corresponding reference simulation (S For simulations using either the composite map or the individual cultivation maps, the differences in the outputs between the reference and bioenergy crop scenarios were analyzed to determine the influence of each bioenergy crop cultivation. Here, we mainly focus on the spatial patterns of bioenergy cultivation impacts, and thus we used the outputs averaged over the last 10 years of each scenario, covering two 5-year rotation cycles set in the model for woody bioenergy crop harvest. Decomposing ΔT a into different factors. Changes of air temperature at 2-m height (ΔT a ) in response to bioenergy crop cultivation were decomposed into components that represent key processes and mechanisms of air temperature change following the methods of Li et al. 24 and Zeng et al. 23 . The decomposition is fundamentally based on the land-surface energy balance: where SW dn , SW up , LW dn , and LW up , denote downward shortwave radiation, upward shortwave radiation, downward longwave radiation, and upward longwave radiation, respectively. LE, H, and G represent the latent, sensible, and ground heat fluxes, respectively. Net shortwave radiation at the surface (SW dn −SW up ) can be expressed as a combination of surface albedo (α) and downward shortwave radiation (SW dn ): Net longwave radiation at the surface (LW dn −LW up ) can be expressed as a function of air temperature (T a ), surface temperature (T s ), and the longwave radiative emissivity of air (ε a ) according to the Stephan-Boltzmann law: where ε s is surface longwave radiative emissivity (prescribed as 1.0 in the coupled model) and σ is the Stephan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ). Sensible heat flux (H) can be given by: where r a , ρ, and C P represent aerodynamic resistance, air density, and the specific heat capacity of air, respectively. Based on Eqs. (4)-(6), Eq. (3) can be re-written as: ε s σT 4 s þ ρC P r a T s ¼ SW dn 1 À α ð ÞÀLE À G þ ε s σε a T 4 a þ ρC P r a T a ð7Þ Differentiation with respect to T s gives: where f s and f a are the energy redistribution factors for surface temperature and air temperature: f a ¼ 4ε s σε a T 3 a þ ρC P r a ð10Þ ΔT a can be partitioned into changes in the local surface energy balance (ΔT local a ) and changes in atmospheric circulation (ΔT cir a ) 23,24 : Similarly, 4T s can be attributed to ΔT local s and ΔT cir s : The first six terms on the right hand side of Eq. (8) Using Eqs. (8), (11), (13), and (14), changes of air temperature change can be decomposed as: where ÀSW dn Δα, 1 À α ð ÞΔSW dn , ÀΔLE, ÀΔG, ρÁC p r 2 a T s À T a À Á Δr a , and ε s σT 4 a Δε a represent air temperature changes induced by changes in albedo, shortwave downward radiation, latent heat flux (mainly evapotranspiration), ground heat flux, aerodynamic resistance, and air emissivity, respectively. Similar to Eq. (13), these six terms constitute the effects of the altered local surface energy balance: Temperature change induced by altered atmospheric circulation (ΔT cir a ) is calculated as the residual from Eq. (11). The altered atmospheric circulation would affect air temperature far away from the regions with land cover change (e.g., the BECCS regions in this study). The relative contribution of ΔT cir a to ΔT a deduced in our study is comparable to previous estimates 23

Data availability
The data which support the main findings of this study are provided as the Source data files of this paper, which is also available via https://doi.org/10.5281/ zenodo.5712916. Source data are provided with this paper.

Code availability
The source code of ORCHIDEE-MICT-BIOENERGY is available at https://doi.org/10.14768/ 02v2-z742, and more details can be found at https://forge.ipsl.jussieu.fr/orchidee/wiki/ GroupActivities/CodeAvalaibilityPublication/ORCHIDEE_MICT_BIOENERGY_r7298. ORCHIDEE-MICT-BIOENERGY is governed by the CeCILL licence under French law and abides by the rules of the distribution of free software. One can use, modify, and/or redistribute the software under the terms of the CeCILL licence as circulated by CEA, CNRS, and INRIA at the following URL: http://www.cecill.info (last access: December 4, 2021).