# Global cooling induced by biophysical effects of bioenergy crop cultivation

Jingmeng Wang, Wei Li, Philippe Ciais, Laurent Z X Li, Jinfeng Chang, Daniel Goll, Thomas Gasser, Xiaomeng Huang, Narayanappa Devaraju, Olivier Boucher

► To cite this version:

Jingmeng Wang, Wei Li, Philippe Ciais, Laurent Z X Li, Jinfeng Chang, et al.. Global cooling induced by biophysical effects of bioenergy crop cultivation. Nature Communications, 2021, 12 (1), 10.1038/s41467-021-27520-0. hal-03501955

HAL Id: hal-03501955

https://hal.science/hal-03501955

Submitted on 24 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-021-27520-0

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Global cooling induced by biophysical effects of bioenergy crop cultivation

Jingmeng Wang<sup>1</sup>, Wei Li<sup>1,2✉</sup>, Philippe Ciais<sup>3</sup>, Laurent Z. X. Li<sup>4</sup>, Jinfeng Chang<sup>5</sup>, Daniel Goll<sup>6</sup>, Thomas Gasser<sup>7</sup>, Xiaomeng Huang<sup>1</sup>, Narayanappa Devaraju<sup>3</sup> & Olivier Boucher<sup>8</sup>

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.

<footnote n="1">¹Department of Earth System Science, Ministry of Education Key Laboratory for Earth System Modeling, Institute for Global Change Studies, Tsinghua University, Beijing 100084, China. ²Institute for Carbon Neutrality, Tsinghua University, Beijing 100084, China. ³Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France. ⁴Laboratoire de Météorologie Dynamique, Centre National de la Recherche Scientifique, Sorbonne Université, Ecole Normale Supérieure, Ecole Polytechnique, 75252 Paris, France. ⁵College of Environmental and Resource Sciences, Zhejiang University, Hangzhou, China. ⁶Université Paris Saclay, CEA-CNRS-UVSQ, LSCE/IPSL, Gif sur Yvette, France. ⁷International Institute for Applied Systems Analysis (IIASA), Laxenburg, Austria. ⁸Institut Pierre-Simon Laplace, Sorbonne Université/IPSL, Paris, France.</footnote>

<footnote>✉email: wli2019@tsinghua.edu.cn</footnote>
Large-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₂ from the atmosphere¹⁻³. Most mitigation scenarios which limit the global temperature increase below 2 or 1.5 °C in 2100 assume widespread deployment of this NET²,⁴. In the scenarios where BECCS was the only land-based negative emission option⁵,⁶, 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⁷, available cultivation lands⁸, land use change emissions⁵,⁶, and supply of water⁹,¹⁰ and nutrients¹¹⁻¹³. 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⁸,¹⁴,¹⁵, but large-scale use of current agricultural systems for bioenergy crops was also proposed¹²,¹⁶,¹⁷. 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¹⁸, but the biophysical interactions between large-scale 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¹⁹⁻²³. In turn, these changes modify air temperature (Tₐ) locally and regionally, in the latter case through atmospheric circulation changes and teleconnection mechanisms²⁴. 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²⁵). IPSL-CM deploys a land surface scheme with an explicit representation of the main lignocellulosic bioenergy crops (ORCHIDEE-MICT-BIOENERGY²⁶) as four plant functional types (PFTs): eucalypt, poplar & willow (sharing the same PFT because of their similar ecological processes and parameters), miscanthus, and switchgrass²⁶. 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²⁶,²⁷. Model performance was evaluated against those field measured yields²⁷ and observation-based 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–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ₐ) 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 composite map (Supplementary Fig. 10f) for bioenergy crop deployment around the world was generated from three marginal land datasets⁸,¹⁴ and two bioenergy cultivation maps from two different IAMs¹⁷ (see Supplementary Fig. 10a–e and Supplementary Methods 3). The

global total bioenergy crop cultivation area in the composite cultivation map is 466 Mha, distributed between 38°S to 60°N (the BECCS regions in Fig. 1a). Based on the composite map, we performed five 50-year simulations with the IPSL-CM model to simulate the biophysical effects of different bioenergy crop cultivation scenarios (more details in Supplementary Methods 4.1 and Supplementary Figs. 10–18): a reference scenario (Sref) where grid cells in the BECCS regions are fully occupied by food crops (because current marginal land areas are low yield croplands, or recently abandoned arable land, see Supplementary Methods 3 and 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 Sref. 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ₐ = -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ₐ 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ₐ range from 0.38 to 0.42 °C across different scenarios, much higher than the global mean ΔTₐ. Strong interannual variability (IAV) was also found in the last 20 years of the global mean ΔTₐ (IAV = 0.16–0.19 °C, see details in Supplementary Methods 4.3, Supplementary Figs. 19 and 20).

Spatial patterns of ΔTₐ differ among the four idealized bioenergy crop scenarios using the composite cultivation map (Fig. 1c-f). Tₐ 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ₐ 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ₐ 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ₐ 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).
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Bioenergy crop cultivation map and air temperature change over land. **a** The composite bioenergy crop cultivation map where red grid cells represent the bioenergy crop cultivation regions (i.e., BECCS regions) and the main vegetation type in each grid cell is shown outside the BECCS regions. **c-f** Air temperature changes ($\Delta T_a$) induced by cultivation of each bioenergy crop type based on the composite cultivation map with their zonal averages shown in (b). Black shading in c-f indicates grid cells with significant ($p < 0.1$) change in air temperature.</figcaption></figure>

**Contribution of energy fluxes to $\Delta T_a$.** We decomposed the $\Delta 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:

$$ \Delta T_a = \Delta T_a^{\text{local}} + \Delta T_a^{\text{cir}} \quad (1) $$

where $\Delta T_a^{\text{local}}$ 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). $\Delta T_a^{\text{cir}}$ represents the impact of atmospheric circulation changes induced by changing the coverage of bioenergy crops (see details in “Methods” and Supplementary Methods 2.4). $\Delta T_a^{\text{local}}$ can be determined from refs. 23,24:

$$ \Delta T_a^{\text{local}} = \frac{1}{f_a} \left[ -\text{SW}_{\text{dn}}\Delta\alpha + (1-\alpha)\Delta\text{SW}_{\text{dn}} - \Delta\text{LE} - \Delta G + \frac{\rho \cdot C_p}{r_a^2} (T_s - T_a)\Delta r_a + \epsilon_s \sigma T_a^4 \Delta\epsilon_a \right] \quad (2) $$

where $f_a$ is an energy redistribution factor (see “Methods”), and $\Delta SW_{\text{dn}}$, $\Delta\alpha$, $\Delta LE$, $\Delta G$, $\Delta r_a$, and $\Delta\epsilon_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$, $\rho$, $\epsilon_s$, and $\sigma$ 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 \times 10^{-8} \text{ W} \cdot \text{m}^{-2} \cdot \text{K}^{-4}$), respectively.

Globally, for woody bioenergy crop scenarios based on the composite cultivation map (i.e., $S_{euc}$ and $S_{p&w}$), $\Delta T_a^{\text{local}}$ 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 $\Delta T_a^{\text{cir}}$ differs between $S_{euc}$ (warming) and $S_{p&w}$ (cooling) (Fig. 2a). Compared to $S_{p&w}$, the cooling effect of $\Delta T_a^{\text{local}}$ in $S_{euc}$ is largely counteracted by the warming effect from $\Delta T_a^{\text{cir}}$, leading to a smaller net $\Delta T_a$. In contrast, for herbaceous bioenergy crop scenarios, $\Delta T_a^{\text{local}}$ is a net warming effect, while $\Delta T_a^{\text{cir}}$ is a cooling effect (Fig. 2a). For both $S_{mis}$ and $S_{swi}$, the cooling effect of $\Delta T_a^{\text{cir}}$ outweighs the warming effect of $\Delta T_a^{\text{local}}$ causing a global net cooling effect (global mean $\Delta T_a < 0$). The warming effect of $\Delta T_a^{\text{local}}$ for the herbaceous bioenergy crop scenarios is mainly contributed by increased aerodynamic resistance, decreased latent heat and increased downward short-wave 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_a^{\text{local}}$ (Fig. 2b and Supplementary Figs. 14–17), whereas $S_{swi}$ decreases evapotranspiration which offsets the cooling effects from other terms in Eq.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Contributions of different components to air temperature changes. Three panels represent the results at the global scale (a), in the BECCS regions (b), 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 $\Delta T_a^{\text{cir}}$, G, $\varepsilon_a$, $r_a$, LE, $SW_{dn}$, and $\alpha$ 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.</figcaption></figure>

(2) (e.g., a higher albedo), resulting in an increased $\Delta T_a^{\text{local}}$ (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, $\Delta T_a^{\text{local}}$ shows warming effects for all the idealized bioenergy crop scenarios based on the composite cultivation map (Fig. 2c), while $\Delta T_a^{\text{cir}}$ 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_{\text{ref}}$), atmospheric circulation and teleconnections alter local radiative forcing indirectly through exchanges of heat between different regions (Supplementary Figs. 14–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**<sub>a</sub> changes at continental scale. Based on the composite bioenergy crop cultivation map, the effects of local surface energy balance changes in air temperature ($\Delta T_a^{\text{local}}$) range from warming in Oceania, which has the smallest bioenergy crop cultivation area, to
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Bioenergy crop cultivation area and mean temperature change in each continent. Air temperature changes induced by altered local surface energy balance (a) and by atmospheric circulation (b) and the net air temperature change (c) are the average values over the whole continent. Sizes of rectangles represent relative bioenergy crop cultivation areas in different continents in the composite cultivation map, and color gradients indicate temperature change. Numbers show the area of bioenergy crop cultivation in each continent (left, M ha) in the composite cultivation map and its proportion of the total continental area (right, %).</figcaption></figure>

cooling in Europe, which has the largest area of bioenergy crop cultivation (Fig. 3a). In fact, $\Delta T_a^{local}$ 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 ($\Delta T_a^{cir}$, Fig. 3b) that redistribute energy spatially show a more subtle pattern (Fig. 3b), either enhancing or weakening the $\Delta T_a^{local}$ signal. As a consequence, mean $\Delta T_a$ at the continental scale depends on both $\Delta T_a^{local}$ and $\Delta T_a^{cir}$ 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_a^{local}$ in Europe ($\Delta T_a^{local} = -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

$\Delta T_a$ is also contributed by $\Delta T_a^{cir}$ (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 $\Delta T_a^{local}$ is counteracted by a warming effect of $\Delta T_a^{cir}$, 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 $\Delta 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 $\Delta T_a^{local}$ and $\Delta T_a^{cir}$ 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 ($\Delta T_a > 0$) contributed by positive $\Delta T_a^{cir}$ (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), but the former magnitudes are smaller (Supplementary Figs. 26 and 29).

**Sensitivity tests of additional scenarios.** To quantify the sensitivity of biophysical climate effects to the choice of a bioenergy crop cultivation map, we performed six additional experimental simulations (Supplementary Methods 4.2), comprising three individual cultivation maps (two IAM maps¹⁷ from IMAGE and MAgPIE, and one marginal land map from Campbell et al.¹⁴, namely Campbell-high, 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 $\Delta T_a^{\text{local}}$ and most energy components for eucalypt are robust among all cultivation maps, but the magnitudes vary, especially for $\Delta T_a^{\text{cir}}$ (Fig. 4). Different from the net cooling effect of $\Delta T_a$ using composite, IMAGE, or MAgPIE map, the positive $\Delta T_a$ in $S_{\text{euc}}^{\text{Cam}}$ (i.e., the simulation using the Campbell-high map for eucalypt) is caused by the large warming effect of $\Delta T_a^{\text{cir}}$, which offsets the cooling effect of $\Delta T_a^{\text{local}}$. By contrast, the changes of biophysical effects of switchgrass are different among simulations with different cultivation maps. $\Delta T_a$ shows a cooling effect in $S_{\text{swi}}$, but a warming effect in $S_{\text{swi}}^{\text{IMA}}$ and $S_{\text{swi}}^{\text{Mag}}$, and a nearly neutral effect in $S_{\text{swi}}^{\text{Cam}}$. Compared to eucalypt, the responses of each component to switchgrass cultivations are relatively small using all cultivation maps, resulting in higher sensitivity of $\Delta T_a$ to cultivation maps. The signs of $\Delta T_a$ for switchgrass are largely determined by $\Delta T_a^{\text{cir}}$.

Similarly, at the continental scale, the signal of $\Delta T_a^{\text{local}}$ 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 $\Delta T_a$ also depend on the $\Delta T_a^{\text{cir}}$ (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 $\Delta T_a^{\text{local}}$ in Europe and $\Delta T_a^{\text{cir}}$ in Oceania and Africa are consistent for switchgrass cultivation across the four maps (one composite and three individual maps). Combining $\Delta T_a^{\text{local}}$ and $\Delta T_a^{\text{cir}}$, however, no consistent $\Delta 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 $\Delta T_a^{\text{local}}$ and $\Delta T_a$ in $S_{\text{euc}}^{\text{IMA}}$ are negative (Fig. 4), indicating the stronger cooling effects of eucalypt plantations than previous forest types. Compared to $S_{\text{euc}}$

<figure><img src="image_7.png" /><figcaption>**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 $\Delta T_a^{\text{cir}}$, G, $\varepsilon_a$, $r_a$, LE, $SW_{dn}$, and $\alpha$ 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.</figcaption></figure>
and $S_{\text{euc}}^{\text{MAg}}$, where eucalypt was cultivated on short-vegetation lands, the cooling effects in $S_{\text{euc}}^{\text{IMA}}$ 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_{\text{swi}}^{\text{IMA}}$, both $\Delta T_{\text{a}}^{\text{local}}$ and $\Delta T_{\text{a}}$ show warming effects ($\Delta T_{\text{a}}^{\text{local}} = 0.01$ °C, $\Delta T_{\text{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 $\Delta T_{\text{a}}^{\text{cir}}$ 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 $\Delta T_{\text{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⁶, 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³⁰,³¹, 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 $\Delta T_{\text{a}}$ and $\Delta T_{\text{a}}^{\text{local}}$ (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_{\text{swi}}^{\text{IMA}}$ 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, $\Delta T_{\text{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⁶, the water demand⁹,¹⁰ and the nutrient cycle¹¹⁻¹³. A full assessment of the nexus of these aspects is thus needed.

## Methods

**Model and simulations.** To simulate the biophysical feedbacks of large-scale bioenergy crop cultivation, ORCHIDEE-MICT-BIOENERGY was coupled with the atmospheric general circulation model (LMDz)³²,³³ of the Institute Pierre Simon Laplace climate model (IPSL-CM²⁵, see Supplementary Methods 1 and 2 for details). The ocean model was excluded from the coupling, and sea surface temperature (SST) was prescribed using climatology data with seasonal variations (from the Atmospheric Model Intercomparison Project (AMIP), http://www-pcmdi.llnl.gov/projects/amip). The spatial resolution of the coupled model was 1.26° latitude × 2.5° longitude.

ORCHIDEE-MICT-BIOENERGY was developed specifically to represent bioenergy crop modeling²⁶ based on ORCHIDEE-MICT³⁴. Two woody and two herbaceous bioenergy crop types were implemented: (1) eucalypts, (2) poplar & willow (treated as the same PFT), (3) miscanthus, and (4) switchgrass.

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²⁶,²⁷. The model is capable of capturing the growth dynamics of bioenergy crops and reproducing the biomass yields in a large global yield observation dataset²⁶,²⁷. The LMDz used here is version 6 of the Laboratoire de Météorologie Dynamique atmospheric general circulation model with zooming capability³²,³³. 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_{\text{ref}}$) and four idealized bioenergy crop scenarios with cultivation of different bioenergy crops ($S_{\text{euc}}$: eucalypt, $S_{\text{p&w}}$: poplar & willow, $S_{\text{mis}}$: miscanthus, $S_{\text{swi}}$: switchgrass, see Supplementary Methods 4.1). In $S_{\text{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_{\text{ref}}^{\text{pre}}$) and six more realistic bioenergy crop scenarios (i.e., 3 individual cultivation maps × 2 representative bioenergy crop types, Supplementary Methods 4.2 and Supplementary 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_{\text{ref}}^{\text{pre}}$) was also run to match the six additional simulations. Different from $S_{\text{ref}}$ with
food crop in the bioenergy crop cultivation grid cells, $S_{ref}^{pre}$ used the present-day observed fractional land covers in each grid cell (Supplementary Methods 4.2).

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 crop 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<sub>a</sub> into different factors.** Changes of air temperature at 2-m height (ΔT<sub>a</sub>) 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.<sup>24</sup> and Zeng et al.<sup>23</sup>. The decomposition is fundamentally based on the land-surface energy balance:

$$SW_{dn} - SW_{up} + LW_{dn} - LW_{up} = LE + H + G \quad (3)$$

where SW<sub>dn</sub>, SW<sub>up</sub>, LW<sub>dn</sub>, and LW<sub>up</sub>, 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<sub>dn</sub>–SW<sub>up</sub>) can be expressed as a combination of surface albedo (α) and downward shortwave radiation (SW<sub>dn</sub>):

$$SW_{dn} - SW_{up} = SW_{dn} (1 - \alpha) \quad (4)$$

Net longwave radiation at the surface (LW<sub>dn</sub>–LW<sub>up</sub>) can be expressed as a function of air temperature (T<sub>a</sub>), surface temperature (T<sub>s</sub>), and the longwave radiative emissivity of air (ε<sub>a</sub>) according to the Stephan–Boltzmann law:

$$LW_{dn} - LW_{up} = \varepsilon_s \sigma (\varepsilon_a T_a^4 - T_s^4) \quad (5)$$

where ε<sub>s</sub> is surface longwave radiative emissivity (prescribed as 1.0 in the coupled model) and σ is the Stephan–Boltzmann constant (5.67 × 10<sup>-8</sup> W·m<sup>-2</sup>·K<sup>-4</sup>).

Sensible heat flux (H) can be given by:

$$H = \rho \cdot C_p \frac{T_s - T_a}{r_a} \quad (6)$$

where r<sub>a</sub>, ρ, and C<sub>p</sub> 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:

$$\varepsilon_s \sigma T_s^4 + \frac{\rho C_p}{r_a} T_s = SW_{dn} (1 - \alpha) - LE - G + \varepsilon_s \sigma \varepsilon_a T_a^4 + \frac{\rho C_p}{r_a} T_a \quad (7)$$

Differentiation with respect to T<sub>s</sub> gives:

$$\Delta T_s = \frac{1}{f_s} \left[ -SW_{dn}\Delta\alpha + (1 - \alpha)\Delta SW_{dn} - \Delta LE - \Delta G + \frac{\rho \cdot C_p}{r_a^2} (T_s - T_a)\Delta r_a + \varepsilon_s \sigma T_a^4 \Delta \varepsilon_a \right] + \frac{f_s}{r_a} \Delta T_a \quad (8)$$

where f<sub>s</sub> and f<sub>a</sub> are the energy redistribution factors for surface temperature and air temperature:

$$f_s = 4\varepsilon_s \sigma T_s^3 + \frac{\rho C_p}{r_a} \quad (9)$$

$$f_a = 4\varepsilon_s \sigma \varepsilon_a T_a^3 + \frac{\rho C_p}{r_a} \quad (10)$$

ΔT<sub>a</sub> can be partitioned into changes in the local surface energy balance (ΔT<sub>a</sub><sup>local</sup>) and changes in atmospheric circulation (ΔT<sub>a</sub><sup>cir</sup>)<sup>23,24</sup>:

$$\Delta T_a = \Delta T_a^{\text{local}} + \Delta T_a^{\text{cir}} \quad (11)$$

Similarly, ΔT<sub>s</sub> can be attributed to ΔT<sub>a</sub><sup>local</sup> and ΔT<sub>a</sub><sup>cir</sup>:

$$\Delta T_s = \Delta T_s^{\text{local}} + \Delta T_s^{\text{cir}} \quad (12)$$

The first six terms on the right hand side of Eq. (8) constitute ΔT<sub>a</sub><sup>local</sup>:

$$\Delta T_s^{\text{local}} = \frac{1}{f_s} \left[ -SW_{dn}\Delta\alpha + (1 - \alpha)\Delta SW_{dn} - \Delta LE - \Delta G + \frac{\rho \cdot C_p}{r_a^2} (T_s - T_a)\Delta r_a + \varepsilon_s \sigma T_a^4 \Delta \varepsilon_a \right] + \Delta T_a^{\text{cir}} \quad (13)$$

ΔT<sub>s</sub><sup>local</sup> is further converted to ΔT<sub>a</sub><sup>local</sup> by Li et al.<sup>24</sup> and Zeng et al.<sup>23</sup>:

$$\Delta T_a^{\text{local}} = \frac{f_s}{f_a} \Delta T_s^{\text{local}} \quad (14)$$

Using Eqs. (8), (11), (13), and (14), changes of air temperature can be decomposed as:

$$\Delta T_s = \frac{1}{f_s} \left[ -SW_{dn}\Delta\alpha + (1 - \alpha)\Delta SW_{dn} - \Delta LE - \Delta G + \frac{\rho \cdot C_p}{r_a^2} (T_s - T_a)\Delta r_a + \varepsilon_s \sigma T_a^4 \Delta \varepsilon_a \right] + \Delta T_a^{\text{cir}} \quad (15)$$

where −SW<sub>dn</sub>Δα, (1 − α)ΔSW<sub>dn</sub>, −ΔLE, −ΔG, ρ·C<sub>p</sub>/r<sub>a</sub><sup>2</sup> (T<sub>s</sub> − T<sub>a</sub>)Δr<sub>a</sub>, and ε<sub>s</sub>σT<sub>a</sub><sup>4</sup>Δε<sub>a</sub> 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:

$$\Delta T_a^{\text{local}} = \frac{1}{f_a} \left[ -SW_{dn}\Delta\alpha + (1 - \alpha)\Delta SW_{dn} - \Delta LE - \Delta G + \frac{\rho \cdot C_p}{r_a^2} (T_s - T_a)\Delta r_a + \varepsilon_s \sigma T_a^4 \Delta \varepsilon_a \right] \quad (16)$$

Temperature change induced by altered atmospheric circulation (ΔT<sub>a</sub><sup>cir</sup>) 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<sub>a</sub><sup>cir</sup> to ΔT<sub>a</sub> deduced in our study is comparable to previous estimates<sup>23,35</sup> (see Supplementary Methods 2.4).

## 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-742, and more details can be found at https://forge.ipsl.jussieu.fr/orchidee/wiki/GroupActivities/CodeAvailabilityPublication/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).

Received: 25 January 2021; Accepted: 23 November 2021;

Published online: 13 December 2021

## References

1. Hepburn, C. et al. The technological and economic prospects for CO₂ utilization and removal. Nature 575, 87–97 (2019).

2. Rogelj, J. et al. Scenarios towards limiting global mean temperature increase below 1.5 °C. Nat. Clim. Change 8, 325–332 (2018).

3. Rosen, J. The carbon harvest. Science 359, 733–737 (2018).

4. Smith, P. et al. Biophysical and economic limits to negative CO₂ emissions. Nat. Clim. Change 6, 42–50 (2016).

5. Harper, A. B. et al. Land-use emissions play a critical role in land-based mitigation for Paris climate targets. Nat. Commun. 9, 2938 (2018).

6. Krause, A. et al. Large uncertainty in carbon uptake potential of land-based climate-change mitigation efforts. Glob. Change Biol. 24, 3025–3038 (2018).

7. Li, W. et al. Mapping the yields of lignocellulosic bioenergy crops from observations at the global scale. Earth Syst. Sci. Data 12, 789–804 (2020).

8. Cai, X., Zhang, X. & Wang, D. Land availability for biofuel production. Environ. Sci. Technol.: Est. 45, 334–339 (2011).

9. Næss, J. S., Cavalett, O. & Cherubini, F. The land-energy-water nexus of global bioenergy potentials from abandoned cropland. Nat. Sustain. 4, 525–536 (2021).

10. Stenzel, F. et al. Irrigation of biomass plantations may globally increase water stress more than climate change. Nat. Commun. 12, 1512 (2021).

11. Bouwman, A. F., Van Grinsven, J. J. & Eickhout, B. Consequences of the cultivation of energy crops for the global nitrogen cycle. Ecol. Appl. 20, 101–109 (2010).

12. Humpenoeder, F. et al. Large-scale bioenergy production: how to resolve sustainability trade-offs? Environ. Res. Lett. 13, 24011 (2018).

13. Li, W. et al. Bioenergy crops for low warming targets require half of the present agricultural fertilizer use. Environ. Sci. Technol. 55, 10654–10661 (2021).

14. Campbell, J. E., Lobell, D. B., Genova, R. C. & Field, C. B. The global potential of bioenergy on abandoned agriculture lands. Environ. Sci. Technol. 42, 5791–5794 (2008).

15. Nijsen, M., Smeets, E., Stehfest, E. & Vuuren, D. P. An evaluation of the global potential of bioenergy production on degraded lands. Gcb. Bioenergy. 4, 130–147 (2012).

16. Perlack, R. D. et al. Biomass as Feedstock for a Bioenergy and Bioproducts Industry: the Technical Feasibility of a Billion-ton Annual Supply (2005).

17. Krause, A., Bayer, A. D., Pugh, T. Å. M., Doelman, J. C. & Arneth, A. Global consequences of afforestation and bioenergy cultivation on ecosystem service indicators. Biogeosciences 14, 4829–4850 (2017).

18. Georgescu, M., Lobell, D. B. & Field, C. B. Direct climate effects of perennial bioenergy crops in the United States. Proc. Natl Acad. Sci. USA 108, 4307–4312 (2011).
19. Alkama, R. & Cescatti, A. Biophysical climate impacts of recent changes in global forest cover. *Science* **351**, 600–604 (2016).

20. Betts, R. A. Offset of the potential carbon sink from boreal forestation by decreases in surface albedo. *Nature* **408**, 187–190 (2000).

21. Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. *Nature* **479**, 384–387 (2011).

22. Luyssaert, S. et al. Land management and land-cover change have impacts of similar magnitude on surface temperature. *Nat. Clim. Change* **4**, 389–393 (2014).

23. Zeng, Z. et al. Climate mitigation from vegetation biophysical feedbacks during the past three decades. *Nat. Clim. Change* **7**, 432–436 (2017).

24. Li, Y., Piao, S., Chen, A., Ciais, P. & Li, L. Z. X. Local and teleconnected temperature effects of afforestation and vegetation greening in China. *Natl Sci. Rev.* **7**, 897–912 (2020).

25. Boucher, O. et al. Presentation and evaluation of the IPSL-CM6A-LR climate model. *J. Adv. Model. Earth Sy.* **12**, e2019MS002010 (2020).

26. Li, W. et al. ORCHIDEE-MICT-BIOENERGY: an attempt to represent the production of lignocellulosic crops for bioenergy in a global vegetation model. *Geosci. Model Dev.* **11**, 2249–2272 (2018).

27. Li, W., Ciais, P., Makowski, D. & Peng, S. A global yield dataset for major lignocellulosic bioenergy crops based on field measurements. *Sci. Data* **5**, 180169 (2018).

28. Abigail, L. S. S., Fung, I. Y. & John, C. H. C. Mid-latitude afforestation shifts general circulation and tropical precipitation. *Proc. Natl Acad. Sci. USA* **109**, 712–716 (2012).

29. Lawrence, D. & Vandecar, K. Effects of tropical deforestation on climate and agriculture. *Nat. Clim. Change* **5**, 27–36 (2015).

30. Zona, D. Long-term effects of permafrost thaw. *Nature* **537**, 625–626 (2016).

31. Schuur, E. A. G. et al. Climate change and the permafrost carbon feedback. *Nature* **520**, 171–179 (2015).

32. Hourdin, F. et al. The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection. *Clim. Dynam.* **27**, 787–813 (2006).

33. Contoux, C., Ramstein, G. & Jost, A. Modelling the mid-Pliocene Warm Period climate with the IPSL coupled model and its atmospheric component LMDZ5A. *Geosci. Model Dev.* **5**, 903–917 (2012).

34. Guimberteau, M. et al. ORCHIDEE-MICT (v8.4.1), a land surface model for the high latitudes: model description and validation. *Geosci. Model Dev.* **11**, 121–163 (2018).

35. Luyssaert, S. et al. Trade-offs in using European forests to meet climate objectives. *Nature* **562**, 259–262 (2018).

## Acknowledgements

This study was funded by the National Natural Science Foundation of China (grant number: 42175169, to W.L.), the National Key R&D Program of China (grant number: 2019YFA0606604, to W.L.), and Tsinghua University Initiative Scientific Research

Program (grant number: 20213080023, to W.L.). P.C. and D.G. acknowledge support from the European Research Council Synergy project SyG-2013-610028 IMBALANCE-P (to P.C. and D.G.) and the ANR CLAND Convergence Institute (to P.C. and D.G.). This study was also supported by the National Key Scientific and Technological Infrastructure project "Earth System Numerical Simulation Facility" (i.e., EarthLab, to X.H.).

## Author contributions

W.L. and P.C. designed the study. W.L., X.H., N.D., and J.W. carried out modeling and analysis. J.W. and W.L. wrote the first draft. J.W., W.L., P.C., L.L., J.C., D.G., T.G., and O.B. contribute to the interpretation of the results. All coauthors contributed to the draft revision.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-27520-0.

**Correspondence** and requests for materials should be addressed to Wei Li.

**Peer review information** *Nature Communications* thanks Joyce Bosmans, Daniel Sanchez and the other anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing,

adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021