Cumulative emissions of anthropogenic carbon dioxide (CO2) have been driving the long-term warming1–3, which has predominantly negative impacts on the physical environment, ecosystem, and humanity4,5. To minimize the potential risks of climate change, the 2015 Paris Agreement aims to keep the temperature rise well below 2 °C and make efforts to limit the increase to 1.5 °C, compared to the pre-industrial level6. In models’ pathways to limit global warming to 1.5 °C, net-zero global CO2 emissions need to be achieved by 2050, and net-negative emission should be achieved to return global warming to 1.5 °C, following a temperature overshoot7. To accomplish this, anthropogenic emissions must be reduced, and carbon dioxide removal (CDR), which removes and permanently stores CO2 from the atmosphere to neutralize anthropogenic emissions and eventually enable net-negative emission, is inevitable even in the case of aggressive mitigation rates7–10.
Despite the increasing attention to CDR in political and economic discussions, there remain many uncertainties related to the effectiveness of CDR to decline temperatures after a peak due to the lack of understanding of the future behavior of the Earth system11–14. The terrestrial biosphere in the Earth system is currently a natural carbon (C) sink, which removes a large fraction of anthropogenic CO2 from the atmosphere15 and can make a valuable contribution to climate change mitigation through afforestation/reforestation16–18. The future terrestrial carbon cycle is a key factor in determining the future levels of CO2 and, hence, climate change. However, there remain many uncertainties in the future behavior of the terrestrial carbon cycle19–21. In particular, the carbon cycle response under net-negative emissions is highly uncertain; thus, there is an urgent need for a better understanding of land C fluxes after reaching net-zero.
So far, studies focusing on the reversibility and hysteresis of land C pool have been undertaken based on idealized CO2 ramp-up and -down forcing experiments using a few Earth system models (ESMs)22–24. Despite their similar experimental designs (the atmospheric CO2 increases at a rate of 1% yr−1 from pre-industrial levels to a quadrupling for 140 years and then decreases at the same rate), the responses of land C stock to CO2 forcing vary among ESMs, probably due to the differences in both climate sensitivities and models’ terrestrial processes. For example, in previous studies, there are differences in the temporal evolution of land C stock, its spatial distribution, the amount of land C stock, and the timing of its peak22–24. A previous study reported that land C stores are largely reversible within the timescale of changing CO2, continue to decrease after CO2 returns to its initial value, and consequently lose ~40 Pg C at the end of the simulation24. However, in other studies22,23, the terrestrial biosphere continues to uptake C immediately after the start of CO2 ramp-down with a significant hysteresis lag due to the inertia of vegetation and as a result stores more C than its initial state.
These different results among ESMs increase the uncertainty in the estimation of land C fluxes after emissions reach net-zero, making it difficult to establish effective climate mitigation strategies, thereby highlighting the necessity of multimodel approaches. Therefore, to comprehensively understand the reversibility and hysteresis of terrestrial C flux and stock in a multimodel context, we analyzed seven ESMs from the Coupled Model Intercomparison Project Phase 6 (CMIP6)25, which performed the Carbon Dioxide Removal Model Intercomparison Project (CDRMIP)12 climate and carbon cycle reversibility (short name: CDR-reversibility) experiment (see Methods and Supplementary Table 1). In this simulation, CO2 concentration is prescribed to increase at 1% yr−1 to the quadrupling and then decreases at the same rate until the pre-industrial CO2 levels, after which the simulation continues for at least 60 years.
Figure 1 shows the temporal evolutions of key variables corresponding to the ramp-up and -down CO2 forcing. The global mean land surface air temperature (SAT) increases by ~6.6 K and peaks at model Year 144, after which it decreases more slowly than in the ramp-up period (Fig. 1a). The land temperature shows a delayed response due to the thermal inertia of the ocean26,27, remaining ~1 K higher than its initial state at the end of CO2 ramp-down and ~0.9 K higher for the next 60 years. The temporal evolution of land precipitation (PRCP) follows the land temperature28,29, showing a delayed peak of about 4 years than the temperature and a greater response on the ramp-down CO2 pathway, likewise the temperature. However, the land PRCP anomaly exhibits a large internal variability, showing an inter-model spread over a wide range.
In addition to these changes in the climate system, an increase in CO2 concentration significantly affects the land carbon fluxes by regulating terrestrial ecosystem processes. At a global scale, the net primary production (NPP), net carbon uptake by vegetation, including autotrophic respiration costs, linearly increases and then returns to its initial value, showing an almost reversible response (Fig. 1b). This is because the global NPP change is mainly caused by the CO2 fertilization effect, i.e., an increase in the rate of photosynthesis resulting from increased CO2 levels30–32. However, heterotrophic respiration (Rh) exhibits a lagged response to CO2 forcing with a delayed peak of about 7 years compared to the CO2 peak and remains higher than its pre-industrial level at the end of the CO2 ramp-down. As the microbial activity is enhanced under warmer and wetter conditions33–35, the hysteresis of the physical climate system causes this delayed Rh response.
The net atmosphere-to-land carbon flux, net biome productivity (NBP), is mainly controlled by the imbalance between NPP and Rh. The NBP exhibits a large hysteresis behavior. The NBP is positive during the CO2 ramp-up period, demonstrating a well-known role of the land as a carbon sink (NPP > Rh)15. Though the amount of C uptake by the land system rapidly increases in the initial stage of the ramp-up period, it soon becomes almost steady despite the gradual increase of the CO2 concentration. Interestingly, the NBP is still positive after the CO2 concentration begins to decrease (for about 23 years in Fig. 1c). That is, the terrestrial biosphere uptakes C for decades after the CO2 peak (total amount: ~63 Gt C). Although the CO2 is prescribed in the present modeling experiments, this result suggests that the terrestrial ecosystem will further contribute to the reduction of CO2 concentration for decades after the net-zero emission is achieved, which would lessen the reliance on the net negative emission.
After that, the land is converted into a C source (NPP < Rh) as the delayed climate system response triggers the hysteresis lag of Rh, NPP immediately returns to its original state. Unlike the ramp-up period, NBP gradually decreases during the CO2 ramp-down period and thereby shows the maximum negative value at the end of the ramp-down period. This result implies that after the peak of atmospheric CO2 concentration, the lower the CO2 concentration, the more intensive the C reduction needed to reduce the atmospheric CO2 concentration due to the C emitted from land. During the restoring period, NBP starts to return to its initial state since NPP does not reduce any more, but the Rh still decreases. The evolution of the NBP response to the ramp-up and -down CO2 forcing leads to hysteresis in the total land C stock, which is a temporal accumulation of NBP. The total land C stock anomaly continues to increase immediately the CO2 ramp-up starts, and it exhibits a significant hysteresis lag (Fig. 1d). Consequently, the land retains more C than its pre-industrial level until the end of the simulation, including a restoring period, indicating its positive role in mitigating anthropogenic climate changes.
Figure 2 shows the latitudinal dependency of the temporal evolution of the physical climate system and terrestrial carbon cycle. The land temperature shows a delayed response to CO2 forcing with a similar scale of delay time in the entire latitudes, remaining warmer than the pre-industrial period for the entire simulation period. Precipitation in the northern mid-high latitudes shows a similar evolution with the global mean precipitation (Fig. 1a), but in 60°S–20°N, the precipitation anomaly is heterogeneous. In the northern mid-high latitudes, the hysteresis behavior of the climate system results in warmer and wetter conditions on the ramp-down pathway even when the CO2 concentration is the same. This change in the climate condition enhances photosynthesis and lengthens the growing season in northern high latitudes where cold temperature limits vegetation growth36,37. Thus, the response of the leaf area index (LAI), an indicator of vegetation growth, to CO2 forcing is delayed as the latitude is higher. On the other hand, in the tropics where the temperature is close to the optimal temperature for photosynthesis38, LAI is almost reversible within the timescale of changing CO2 because the fertilization impact of CO2 concentration is dominant.
This latitudinal dependency of the terrestrial biosphere is more evident in the evolution of NBP (Fig. 2d). The terrestrial biosphere in high latitudes continues to absorb C for decades after atmospheric CO2 decreases due to the warmer and wetter conditions favorable for sub-Arctic vegetation, and thus, the transition of C sink to the source is more delayed than in tropics. That is, the higher the latitude, the slower the transition of NBP. Accordingly, the vegetation carbon (cVeg) is almost reversible in the tropics at the end of the simulation, whereas the vegetation in mid-high latitudes retains more C than its initial state due to its longer reversibility timescale (Fig. 2e). The C stored in the litter–soil system (the sum of cLitter and cSoil) also exhibits the latitudinal-dependent hysteresis similar to cVeg with a greater time lag to CO2 forcing in the sub-Arctic region (50–70°N) (Fig. 2f). However, this response of cLitter and cSoil to CO2 forcing is generally more delayed than cVeg because of the flow of C from plant biomass to litter and from decaying litter to the soil. The increase in the C stored in the litter–soil system is greatest in high latitudes due to slow decomposition in cold environments, but the largest change of cVeg is in the tropics where the forest density is the highest. In addition, interestingly, cLitter and cSoil are slightly negative in the pan-Arctic region above 60°N, unlike other regions, in the ramp-down period, which could be related to the permafrost process.
It is noted that the multimodel ensemble (MME) mean pattern of land C pool differs from that reported in previous studies using a single model22,24, and this pattern varies among ESMs (Supplementary Figs. 1,2). For example, the amount of change in the C stock and its spatial distribution, such as the maximum C uptake region, are different. Moreover, the peak of the land C stock and the timescale of the lagged response to CO2 forcing vary among models. This may be due to the different soil C representation and inclusion of dynamic vegetation, which simulates shifts in the vegetation as a response to climate change, in only a few ESMs 24. Interestingly, two ESMs coupled with dynamic vegetation (GFDL-ESM4 and UKESM1-0-LL) show a larger hysteresis lag in high latitudes and a slight expansion to northwards (Supplementary Table 2 and Fig. 1). It is also interesting that the C stored in the litter–soil system is negative in the pan-Arctic region in only two ESMs (CESM2 and NorESM2-LM), which include a representation of permafrost processes. Besides, the various ecosystem processes, as represented in ESMs, and the differences in the simulated climate change may be responsible for the inter-model diversity, which causes the poor understanding of the Earth system.
To further understand the regional characteristics of the land C system’s hysteresis behavior and its inter-model diversity, we investigated the spatial pattern of the total land C stock. Figures 3a and b show the total land C stock difference between on the CO2 ramp-down and -up pathways at 2×CO2 and 1×CO2, respectively. Though the CO2 concentration is the same, the total land C stock on the CO2 ramp-down phase distinctly increases, especially in boreal forests, the Maritime Continents, and East Asia. Particularly, boreal forests can store C for a long time owing to their high hysteresis lag (Fig. 3b). However, the differences in the land C stock in Amazon (Supplementary Note 1) and the permafrost regions, especially within the Arctic Circle, are not statistically significant due to the diverse response among ESMs (Figs. 3a,b and Supplementary Figs. 3,4). The inter-model diversity of difference in total land C stock estimated by the coefficient of variation (the standard deviation of the spread divided by the mean) is highest in the Arctic region at both 2×CO2 and 1 × CO2 (Figs. 3c and d), indicating the greatest relative variability in high latitudes. This is because there are two exceptional models (CESM2 and NorESM2-LM), which simulate lower land C stock in the ramp-down period than in the ramp-up period, completely different from the other models (Supplementary Fig. 5).
Therefore, we conducted a more detailed analysis to understand the contrasting terrestrial C stock response to CO2 forcing in the permafrost region located above 60 °N (Fig. 4). CESM2 and NorESM2-LM showed a negative cSoil anomaly at the end of changing CO2 (model Year 280), whereas the other models showed a positive cSoil anomaly (Fig. 4a). The negative cSoil anomaly shown by the two models is attributed to the greater Rh anomalies than NPP anomalies in the entire experimental period. Interestingly, these two models (~1810 Gt C) show about 13 times greater climatology of cSoil than the other ESMs (~140 Gt C). Notably, the two ESMs, coupled with Community Land Model 5, include the representation of the permafrost C pool and its physics (Supplementary Table 2). The CLM has improved its components, such as snow density and insulation, soil layer resolution, and soil hydrology, in the cold region to simulate more realistic permafrost distribution, active layer dynamics, and hydrology39. In particular, the introduction of vertically resolved soil biogeochemistry enables the model to generate large C stocks in the northern high‐latitude permafrost domain, as observed (~1700 Pg C)39–41. As CLM5 includes the simulation of C decomposition in the permafrost zone containing large soil organic C stock, Rh exceeds NPP only in CESM2 and NorESM2-LM, thereby simulating a negative cSoil anomaly. However, soil C stock in the other models is remarkably low compared to the observed value due to the absence of permafrost and related processes, so that they possibly underestimate the soil respiration.
As a result, CESM2 and NorESM2-LM (group A) simulate C release to the atmosphere in the permafrost region, unlike the other models (group B) in the later parts of the CO2 ramp-up period because of the sharp increase in Rh (Fig. 4c and Supplementary Fig. 5). Afterward, during the ramp-down and restoring period, the models in group A simulate C emission continuously due to higher Rh than NPP resulting from activated microbial C decomposition in permafrost under a warmer condition. This result indicates that the Arctic terrestrial could be a C source not a C sink in the course of the ramp-up and -down forcing, implying its accelerating role in global warming. However, the other models show a transition from a C sink to a C source similar to the global mean temporal evolution of NBP but much slower.
Consequently, group B models simulate positive total land C stock over the entire experimental period with a hysteresis lag to the CO2 forcing. In group B models, the terrestrial biosphere serves as a C sink, storing more C (~33 Gt C) than its initial state at the end of the simulation. However, the land C stock anomaly in group A models, including the permafrost representation, show quite different behavior; the total land C stock slightly increases in the early phase of the ramp-up period but gradually decreases afterward. In particular, the total land C stock anomaly becomes negative in model Year 140 when the CO2 concentration is the maximum, and then, it gradually decreases through the ramp-down and restoring periods. Thus, the terrestrial biosphere contains much less C at the end of the restoring period due to the ~38 Gt C loss compared to that in the pre-industrial period, despite the same CO2 concentration, indicating an evident irreversible response to CO2 forcing, which worsens global warming and the reversibility of the climate system.
In this study, we investigated the reversibility and hysteresis behavior of land C fluxes and stock, focusing on carbon cycle responses under negative emissions, based on idealized CO2 ramp-up and -down multimodel simulations. We found that the total land C stock continues to increase for decades, even after CO2 removal starts, and keeps considerably higher levels in the ramp-down period compared to that at the same CO2 level in the ramp-up period, indicating the strong hysteresis behavior of land C cycle. Furthermore, land C stock remains high even after the CO2 concentration returns to its initial level. These results imply the potential mitigation of climate change and the possibility of a decrease in the burden on negative emissions due to the hysteresis behavior of land C flux. Such hysteresis behavior is latitudinally dependent, and the timescale of reversibility is much longer in high-latitude regions because of the inertia of vegetation dynamics under favorable climate conditions. At a regional scale, boreal forests can store C for a longer time than any other region under net-negative emissions. These spatiotemporal characteristics can be considered for establishing an effective strategy for natural climate solutions, such as forest management.
However, we also found that the hysteresis behavior of the terrestrial carbon cycle strongly depends on the representation of the permafrost process. With more realistic permafrost C stock, we found that there would be a strong irreversible process of the terrestrial carbon cycle in the permafrost region, which would hinder global warming mitigation. Notably, this contribution of permafrost has often been omitted in both the mitigation pathway and estimation of the remaining C budgets7,11,14. Our results demonstrate the need to quantify future permafrost C emission and integrate it into the estimation of the remaining C budget. In addition, the considerable permafrost C loss after achieving net-zero emission should be considered in climate policy discussions and decisions.
Furthermore, more careful quantification is needed because the cumulative C emission may be underestimated as no land surface model considers abrupt thaw, the rapid degradation of ice-rich permafrost by thaw lake formation, and erosion42,43. In addition, due to the characteristic of path-dependent permafrost C release44 and the limitation of the idealized experiment, it is necessary to re-estimate C emission from the permafrost region by applying more realistic negative emission scenarios. Due to these current limitations in ESMs and the design of experiments associated with permafrost processes, unexpected additional future C emissions in the pan-Arctic region might occur, which increases the uncertainty about future levels of CO2 and climate change.