Hysteresis of terrestrial carbon cycle to CO 2 ramp-up and -down forcing

To prevent excessive global warming, we have faced a situation to reduce net carbon dioxide (CO 2 ) emissions. However, the behavior of Earth’s terrestrial biosphere under negative emissions is highly uncertain. Herein, we show strong hysteresis in the terrestrial carbon cycle in response to CO 2 ramp-up and -down forcing. Owing to the strong hysteresis lag, the terrestrial biosphere stores more carbon at the end of simulations than at its initial state, lessening the burden on net-negative emissions. This hysteresis is latitudinally dependent, showing a longer timescale of reversibility in high latitudes. Particularly, carbon in boreal forests can be stored for a long time. However, the hysteresis of the carbon cycle in the pan-Arctic region depends on the presence of permafrost processes. That is, unexpected irreversible carbon emissions may occur in permafrost even after achieving net-zero emissions, indicating the importance of permafrost processes, which is highly uncertain based on our current knowledge.

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 preindustrial level 6 . In models' pathways to limit global warming to 1.5 °C, net-zero global CO 2 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 overshoot 7 . To accomplish this, anthropogenic emissions must be reduced, and carbon dioxide removal (CDR), which removes and permanently stores CO 2 from the atmosphere to neutralize anthropogenic emissions and eventually enable net-negative emission, is inevitable even in the case of aggressive mitigation rates [7][8][9][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 system [11][12][13][14] . The terrestrial biosphere in the Earth system is currently a natural carbon (C) sink, which removes a large fraction of anthropogenic CO 2 from the atmosphere 15 and can make a valuable contribution to climate change mitigation through afforestation/reforestation [16][17][18] . The future terrestrial carbon cycle is a key factor in determining the future levels of CO 2 and, hence, climate change. However, there remain many uncertainties in the future behavior of the terrestrial carbon cycle [19][20][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 uxes after reaching net-zero.
So far, studies focusing on the reversibility and hysteresis of land C pool have been undertaken based on idealized CO 2 ramp-up and -down forcing experiments using a few Earth system models (ESMs) [22][23][24] . Despite their similar experimental designs (the atmospheric CO 2 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 CO 2 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 peak [22][23][24] . A previous study reported that land C stores are largely reversible within the timescale of changing CO 2 , continue to decrease after CO 2 returns to its initial value, and consequently lose ~40 Pg C at the end of the simulation 24 . However, in other studies 22,23 , the terrestrial biosphere continues to uptake C immediately after the start of CO 2 ramp-down with a signi cant 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 uxes after emissions reach net-zero, making it di cult to establish effective climate mitigation strategies, thereby highlighting the necessity of multimodel approaches. Therefore, to comprehensively understand the reversibility and hysteresis of terrestrial C ux and stock in a multimodel context, we analyzed seven ESMs from the Coupled Model Intercomparison Project Phase 6 (CMIP6) 25  concentration is prescribed to increase at 1% yr −1 to the quadrupling and then decreases at the same rate until the pre-industrial CO 2 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 CO 2 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 ocean 26,27 , remaining ~1 K higher than its initial state at the end of CO 2 ramp-down and ~0.9 K higher for the next 60 years. The temporal evolution of land precipitation (PRCP) follows the land temperature 28,29 , showing a delayed peak of about 4 years than the temperature and a greater response on the ramp-down CO 2 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 CO 2 concentration signi cantly affects the land carbon uxes 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 CO 2 fertilization effect, i.e., an increase in the rate of photosynthesis resulting from increased CO 2 levels [30][31][32] . However, heterotrophic respiration (Rh) exhibits a lagged response to CO 2 forcing with a delayed peak of about 7 years compared to the CO 2 peak and remains higher than its pre-industrial level at the end of the CO 2 ramp-down. As the microbial activity is enhanced under warmer and wetter conditions [33][34][35] , the hysteresis of the physical climate system causes this delayed Rh response.
The net atmosphere-to-land carbon ux, 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 CO 2 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 CO 2 concentration.
Interestingly, the NBP is still positive after the CO 2 concentration begins to decrease (for about 23 years in Fig. 1c). That is, the terrestrial biosphere uptakes C for decades after the CO 2 peak (total amount: ~63 Gt C). Although the CO 2 is prescribed in the present modeling experiments, this result suggests that the terrestrial ecosystem will further contribute to the reduction of CO 2 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 CO 2 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 CO 2 concentration, the lower the CO 2 concentration, the more intensive the C reduction needed to reduce the atmospheric CO 2 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 CO 2 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 CO 2 ramp-up starts, and it exhibits a signi cant 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 CO 2 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 CO 2 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 growth 36,37 . Thus, the response of the leaf area index (LAI), an indicator of vegetation growth, to CO 2 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 photosynthesis 38 , LAI is almost reversible within the timescale of changing CO 2 because the fertilization impact of CO 2 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 CO 2 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 CO 2 forcing in the sub-Arctic region (50-70°N) (Fig. 2f). However, this response of cLitter and cSoil to CO 2 forcing is generally more delayed than cVeg because of the ow 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 model 22,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 CO 2 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 CO 2 ramp-down and -up pathways at 2×CO 2 and 1×CO 2, respectively. Though the CO 2 concentration is the same, the total land C stock on the CO 2 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 signi cant 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 coe cient of variation (the standard deviation of the spread divided by the mean) is highest in the Arctic region at both 2×CO 2 and 1 × CO 2 (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 CO 2 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 CO 2 (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 hydrology 39 . 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][40][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 CO 2 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 CO 2 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 CO 2 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 CO 2 concentration, indicating an evident irreversible response to CO 2 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 uxes and stock, focusing on carbon cycle responses under negative emissions, based on idealized CO 2 ramp-up and -down multimodel simulations. We found that the total land C stock continues to increase for decades, even after CO 2 removal starts, and keeps considerably higher levels in the ramp-down period compared to that at the same CO 2 level in the ramp-up period, indicating the strong hysteresis behavior of land C cycle.
Furthermore, land C stock remains high even after the CO 2 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 ux. 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 budgets 7,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 quanti cation 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 erosion 42,43 . In addition, due to the characteristic of pathdependent permafrost C release 44 and the limitation of the idealized experiment, it is necessary to reestimate 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 CO 2 and climate change.

CMIP6 data and experimental design
The CDR-reversibility experiment from the CMIP6 CDRMIP 12 was analyzed to investigate the carbon cycle response to large-scale CO 2 removal. This experiment was branched from the 1pctCO2 experiment, in which the CO 2 level increases at a rate of 1% yr −1 from pre-industrial levels to quadrupling for 140 years, from the CMIP6 Diagnostic, Evaluation, and Characterization of Klima (DECK) 25 . Then, a 1% yr −1 removal of CO 2 from the atmosphere is prescribed for 140 years until the pre-industrial CO 2 level is reached and then held for as long as possible (minimum of 60 years) in the 1pctCO2-cdr simulation (minimum of 200 years). Thus, the total length of the CDR-reversibility experiment employed herein is 340 years. We calculated the anomalies of CDR-reversibility simulations using the pre-industrial prescribed CO 2 control simulation (piControl) from DECK as a baseline.
We used seven ESMs (ACCESS-ESM1-5, CanESM5, CESM2, GFDL-ESM4, MIROC-ES2L, NorESM2-LM, and UKESM1-0-LL), which were coupled with the full carbon cycle and performed CDR-reversibility simulations from the CMIP6 archive (Supplementary Tables 1,2). The MME mean was derived by regridding the outputs from ESMs to a common 1° × 1° grid, then averaging them. GFDL-ESM4 does not provide cLitter and cSoil data, UKESM1-0-LL does not provide cLitter data, and NorESM2-LM does not provide precipitation data. Due to this limitation, GFDL-ESM4 is excluded when calculating the MME mean of total land C stock, and the total land C stock in UKESM1-0-LL is calculated as the sum of cSoil and cVeg excluding cLitter. Herein, the bootstrap method was used to test the statistical signi cance of the difference between the experiments. For MME, seven values were randomly selected from seven ESMs with replacements, and then, their average was computed. By repeating this process 1000 times, the con dence intervals were determined, and only signi cant values were shown to indicate the model agreement.
Diagnosing permafrost in the model The permafrost extent in the model was diagnosed using the temperature at the minimum soil depth (D zza ) where the monthly mean variation of soil temperature within a year is less than 0.1°C 45 . If the temperature at D zza is less than 0°C for 2 years or more, that grid cell is assumed to be permafrost.
However, there are some CMIP6 models in which a soil pro le is not deep enough to identify D zza (Supplementary Table 2). For such models, permafrost is assumed to be present in grid cells where the 2- year mean soil temperature of the deepest soil layer is less than 0°C 46 . The extent of permafrost in the model was diagnosed using piControl simulation, except for GFDL-ESM4M, which does not provide soil layer temperature data. In the multimodel context, if the grid cell is diagnosed as permafrost in four or more of the six ESMs, we de ne those grid cells as permafrost regions. The MME mean permafrost extent is almost similar to the boundaries of permafrost extent (>50% coverage) from ESA Climate Change

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.