Contrasting ecosystem constraints on seasonal terrestrial CO2 and mean surface air temperature causality projections by the end of the 21st century

Two centuries of studies have demonstrated the importance of understanding the interaction between air temperature and carbon dioxide (CO2) emissions, which can impact the climate system and human life in various ways, and across different timescales. While historical interactions have been consistently studied, the nature of future interactions and the impacts of confounding factors still require more investigation in keeping with the continuous updates of climate projections to the end of the 21st century. Phase 6 of the Coupled Model Intercomparison Project (CMIP6), like its earlier projects, provides ScenarioMIP multi-model projections to assess the climate under different radiative forcings ranging from a low-end (SSP1–2.6) to a high-end (SSP5–8.5) pathway. In this study, we analyze the localized causal structure of CO2, and near-surface mean air temperature (meanT) interaction for four scenarios from three CMIP6 models using a rigorous multivariate information flow (IF) causality, which can separate the cause from the effect within the interaction (CO2–meanT and meanT–CO2) by measuring the rate of IF between parameters. First, we obtain patterns of the CO2 and meanT causal structures over space and time. We found a contrasting emission-based impact of soil moisture (SM) and vegetation (leaf area index (LAI)) changes on the meanT–CO2 causal patterns. That is, SM influenced CO2 sink regions in SSP1–2.6 and source regions in SSP5–8.5, and vice versa found for LAI influences. On the other hand, they function similarly to constrain the future CO2 impact on meanT. These findings are essential for improving long-term predictability where climate models might be limited.


Introduction
The relationship between near-surface air temperature and carbon dioxide emissions (CO 2 ) has remained a topic of significant concern and research for almost two centuries (Arrhenius 1896, Stips et al 2016. Both observed and modeled increases in temperature have been noted to impact CO 2 variability, while increases in CO 2 have also been shown to be a pivotal contributor to the warming climate (Lacis et al 2010). The impacts on various aspects of the ecosystem (van Nes et al 2015, Deryng et al 2016, Demirhan 2020 and socio-economic patterns (Deryng et al 2016, Appiah et al 2018 have been studied from many perspectives using different approaches. Several studies have sought to understand the causality within this interaction using either controlled numerical simulations (for instance, Friedlingstein et al (2001) and Seneviratne et al (2013)) or data-driven approaches (for instance, Attanasio (2012) and Koutsoyiannis and Kundzewicz (2020)), both of which have their strengths and limitations. Based on these approaches, previous studies have explored the multiple temporal characteristics of the interaction (Faes et al 2017) in historical records, ranging from paleoclimate timescales (Stips et al 2016, Barral et al 2017 to annual and 6 month timescales , and even to daily timescales (Kotz et al 2021). These studies have consistently shown that climate-carbon feedback is fundamentally positive. Given the increasing concern regarding the impact of increased CO 2 emissions on the climate and vice versa, assessing how the causal structures would evolve in the future in a changing climate under various warming instances is of utmost importance to improve the predictability of these impacts where climate models are limited.
The primary activity of Phase 6 of the Coupled Model Intercomparison Project (CMIP6)-the Scenario Model Intercomparison Project-has recently implemented multi-model climate projections based on different pathways of future emissions and land use changes (hereafter referred to as CMIP6 emission scenarios) developed with integrated assessment models (O'Neill et al 2016). In the past, proceedings from earlier phases of this project have not only improved our understanding of possible future climate scenarios but have also served as decisionmaking models for characterizing societal risks and related policies. For instance, Samset et al (2020) used the CMIP5 datasets to investigate the impact of individual climate drivers in mitigating projected global surface temperature evolution and found that anthropogenic CO 2 had the most significant potential to influence it. Furthermore, a more recent study used the CMIP5 records to show that rising temperatures could reduce the yields of some major crops (Zhao et al 2017).
Although there is a tight causal relationship between temperature and CO 2 , their impacts on each other are also influenced by some key parameters such as the sea surface temperature in the ocean (Friedlingstein et al 2001), vegetation (Zhu et al 2016), soil respiration (Bond-Lamberty and Thomson 2010) and soil moisture (SM) over land (Green et al 2019, Humphrey et al 2021. This means that changes in these parameters could also lead to changes in the temperature and CO 2 interaction. These suggest the following: what would be the changes in the global causal hotspot locations under the different projection pathways? And how would changes in key influencing factors modulate them in the future? Along these lines, this study attempts to quantify the temperature-CO 2 causality based on the CMIP6 projections under four pathways of socio-economic developments, mitigation of emissions, and pollution control by the end of the 21st century across the globe from an ensemble of three climate models since 1979. Here, we rely on a rigorously formulated information flow (IF)-based causality tool that can quantitatively evaluate the cause-effect relationship in time series (Liang 2014). Furthermore, asymmetry in IF makes it possible to differentiate causal information from mere correlations that do not have causal sufficiency (Liang 2013). This method has been successfully applied to fields such as finance (Liang 2016), neuroscience (Hristopulos et al 2019), climate science (Bai et al 2017, Hagan et al 2019, Docquier et al 2022, and more importantly, historical temperature-CO 2 feedback analysis (Stips et al 2016). Here, we investigate the impact of near-surface atmospheric CO 2 emissions (CO 2 ) on near-surface mean air temperature (meanT) and vice versa (pixel-wise) under four degrees of emission pathways ranging from a sustainable scenario to a business-as-usual scenario chiefly over land. More importantly, we examine how changes in two important land factorsvegetation and SM-also contribute to the future causality between the projected CO 2 and meanT. This allows us to analyze regional ecosystem resilience to changes in the climate-carbon feedback. To capture the impacts of the two factors on the CO 2 and meanT interactions, we focus on the seasonal mean variations since very long timescales would reduce their impacts and very short timescales would be significantly impacted by feedbacks over the land. Thus, we focus on localized short-term impacts.

CMIP6 datasets
Historical records  from the CMIP6 framework and four future scenarios (2015-2100) for pixel-wise meanT from three CMIP6 models (table S1 shows a list of the models used), based on the availability of all four future scenarios of the model CO 2 outputs at the time of the study, are used (O'Neill et al 2016). Here, we use the near-surface emissiondriven CO 2 outputs to include the interactive effects of carbon cycle feedbacks (Friedlingstein et al 2014) necessary for such causal analyses. The four emission pathways (or scenarios) in the projections used here represent four different radiative forcing targets by the end of the 21st century. The selected scenarios here, called shared socio-economic pathways (SSPs), range from SSP1-2.6, the low end of the forcing pathways, through the medium (SSP2-4.5 and SSP3-7.0 respectively) to the high-end, SSP5-8.5, which has the highest emissions among the scenarios. More details are provided in O'Neill et al (2016). Preprocessing steps for the causality computations are presented in the supplementary material (text S1). All the results are based on the ensemble mean of the selected models.

Causality formalism
The causal inference technique used here is from Liang's IF and causality analysis theory rigorously derived from first principles (e.g. Liang 2014Liang , 2016Liang , 2021 and expressed in terms of sample covariances (Liang 2014). Given two time series X 2 (e.g. CO 2 ) and X 1 (e.g. meanT), the causality from the former to the latter proves to be measurable by the time rate of the flow of information from X 2 to X 1 , which is given by (S1).
In reality, a third (confounding) variable could influence the interaction between X 2 and X 1 (e.g. leaf area index (LAI)). Liang (2021) extended the formalism (S1) to multivariate settings. For a system of d number of time series, the IF from X 2 to X 1 influenced by confounding variables in the system is given as: where C j,d1 is the sample covariance between X j and the derived X 1 using the Euler forward differencing, and det C represents the determinant of the sample covariance. It is evident that once d = 2, (1) reduces to (S1). Here, we use the normalized IF of T 2→1|3,4,..,d taken between −100 and 100 (Liang 2021) which allows us to fairly compare the IF results for different degrees of emission scenarios. More details are given in text S2. Following Liang (2014), we use the Fisher information matrix for significance testing since its inverse gives a covariance matrix with a given significance level at a 5% significance level. When T 2→1|3,4,..,d is significant, X 2 is considered to be causal to X 1 .

Results and discussion
3.1. Causal structure between CO 2 and mean air temperature Figure 1 shows the separated cause-effect structure of the interaction, where the two panels show the IF from CO 2 to meanT (CO 2 -meanT) for both historical and future scenarios. We analyze two projection periods: 2015-2060 (figures 1(b)-(e)) and 2061-2100 (figures 1(f)-(i)). Overall, figures 1(a)-(i) show increasing causal strengths and area coverage along the equator with the evolution of time and increasing radiative forcing scenarios. That is, causal strengths are larger toward the end of the century than at the beginning in SSP5-8.5. The positive IFs suggest that changes in CO 2 increase or amplify meanT anomalies in these regions. The projections also show weakly negative IF values in the SSP1-2.6 (figures 1(b) and (f)) which become stronger over the northern hemisphere (NH). This means the CO 2 impact on SSP1-2.6 reverses from amplifying variabilities in meanT to reducing them by the end of the century due to reduced CO 2 emissions that lead to reduced meanT or milder meanT increases. As we move toward the higher end of the scenarios in figures 1(f)-(i), CO 2 increasingly becomes a source of amplifying meanT changes. However, regions in the NH appear unaffected. This latter part will be explored further in the following subsection. Figure 2 demonstrates patterns of sources and sinks for CO 2 due to meanT changes in the different emission pathways and at different century periods. MeanT by itself cannot be a sink or a source for CO 2 ; therefore, these results are more representative of how meanT changes drive sink-source factors of CO 2 . The historical results in figure 2(a) show that, generally, the NH serves mildly as a sink while saturations over the SH lead to more source regions. The strong positive signals found in the tropics, specifically over the Amazon and Congo basins, could be linked to the impact of land-use change in the region. This has led to tree-cover loss (due to deforestation, forest fires and logging), eventually making the region more sensitive to climate warming (Hansen et al 2013). Forest loss generally increases surface albedo, eventually leading to decreased evapotranspiration and LAI (Li et al 2022). As a result, net warming of the region increases, and CO 2 uptake is reduced.
In the projections, the results of the scenarios show that milder pathways (SSP1-2.6 and SSP2-4.5) do lead to increased CO 2 sinks (figures 2(b), (c), (f) and (g)). However, the sources in SSP2-4.5 are anomalously amplified briefly (figure 2(c)) before more sinks appear (figure 2(g)). Friedlingstein et al (2001) showed that after a continuous increase in CO 2 emissions, there would be a point of CO 2 equilibrium or stabilization over both the land and the ocean in the future. Therefore, the change in SSP2-4.5 could be a reflection of that phenomenon. On the other hand, SSP1-2.6 quickly converges to sink scenarios due to reduced emissions (figures 2(b) and (f)). In the high-end scenarios (SSP3-7.0 and SSP5-8.5), continuous emissions contribute to substantial decreases in the sink regions found in the historical and low-end scenario results, which eventually turn into source regions. As expected, significant temperature increases lead to subsequent emissions of CO 2 , especially over the NH, with the strongest impacts found in SSP5-8.5 (figures 2(e) and (i)). Over the SH, we observe that a crucial region like the Amazon basin, which could otherwise function as a sink, turns into a source by the end of the century. The negative IFs found in emission scenarios may be considered as the resilience of global ecosystems to increasing meanT to stabilize CO 2 emissions, thus modulating the positive climate-carbon feedback.
Although meanT has been shown to correlate with CO 2 in previous studies (Kundzewicz et  In the remainder of the paper, we explore these two factors' effects in more detail, especially to explain the potential reasons for the results obtained. To assess the reliability of our results, a similarity analysis (Lo and MacKinlay 1989) is also provided in figures S1 and S2 to show their agreement within the models in both periods. Generally, the models seem to agree in most regions in both periods, which lends confidence to the results obtained in figures 1 and 2. However, care should be taken when interpreting the results over regions of strong disagreement as they may require further verification.

Constraints of LAI and SM on the CO 2 -meanT causality
From the CMIP6 framework, we have chosen LAI as the indicator for vegetation. Here, to assess the influence of third variables on the coupling, we use the normalized multivariate IF causality to estimate the causalities as shown in figures 3(a1)-(e2) for LAI and figures 4(a3)-(e4) for SM. Figure 3 suggests that the evolution of meanT over land depends primarily on CO 2 emissions in this interaction. Nonetheless, there are signs of the influences of LAI and SM. Similar to what we observed in figure 1, CO 2 amplifies meanT variability more in the higher-end scenarios, and by the end of the century most strongly near the equator, especially over the Amazon and the Congo basins. This could be related to the reported turning point of the Amazon rainforest where positive feedback between vegetation and CO 2 will break down at some point in the future due to oversaturation of the net biospheric uptake of CO 2 (Friedlingstein et al 2001) and unfavorable temperatures (Cox et al 2004, Doughty andGoulden 2008). This has been suggested to already be on the brink even in historical records (Huang et al 2019). A similar situation may apply to the Congo basin. In the high-end scenarios of figures 3(b1)-(e3), significant emissions of CO 2 lead to saturations of uptake capacities over these regions that could have otherwise served as sinks due to decreased photosynthesis (Green et al 2020). Figures 3(b) and (c) also show that in the SSP1-2.6, reductions (negative IFs) in CO 2 -meanT appear to be stronger over the NH, although the strength of these reductions decreases with increasing emission scenarios. The faint reductions found in the NH of the high-end scenarios suggest that in the high-end emissions, LAI and SM continue to function as constraints on amplifying meanT variability over land (figures 3(a1)-(e2)) due to impacts on CO 2 and meanT, hence the CO 2 -meanT coupling. In fact, regions where their constraints are strongest in the SSP1-2.6 causality are where we find the weakest or most insignificant IFs in SSP3-7.0 and SSP5-8.5 CO 2 -meanT coupling of the latter part of the century in figures 1(b)-(e) and 3(b2)-(e2). That is, meanT anomalies would be more significant in the high-end scenarios if there were no constraints on the positive climate-carbon feedback. This could represent ecosystem resilience to changes in the carbon-climate feedback that slows down the impact of CO 2 on meanT throughout the century.

Constraints of LAI and SM on the meanT-CO 2 causality
Figures 4(a1)-(e4) explore the influences of LAI and SM changes on the other side of the separatedinteraction meanT-CO 2 causality. Figures 4(a1)-(e2) indicate places where changes in the LAI function influence meanT-CO 2 causality. Negative IFs indicate that LAI contributes to the reductions in the meanT-CO 2 causality, and positive IFs indicate regions where LAI contributes to amplifying the meanT-CO 2 causality. First, a distinct latitudinal spatial pattern can be observed globally from the historical results and projections. Over the NH, LAI amplifies the coupling variability, especially over the Eurasian region, mainly in the low-end scenarios. This is strongest in SSP1-2.6, which intensifies in the second half of the projected period (figures (b1) and (b2)). In figures 4(b1)-(e2), the strong intensification reduces with increasing emission scenarios. Negative IFs appear over the high-latitude regions in SSP5-8.5, although the positive IF regions appear to have moved into the Mediterranean regions, southern USA and eastern parts of Brazil, possibly owing to increased arid conditions, eventually increasing uncertainties in the meanT-CO 2 causality there. Over the SH, the historical results show positive IF due to the influence of LAI on the causality. However, all the projections show negative IFs appearing to get stronger with time and increasing scenarios. This implies that LAI influences reductions in the variability of the causality and consequently reductions in CO 2 variability due to changes in mean Zhu et al (2016) found that climate change resulted in greening over high latitudes. This greening could explain the increases in the negative IFs since greening would also increase photosynthesis and eventually increase CO 2 sinks, as found in the results.
The impact of SM on the coupling, as shown in figures 4(c) and (d), appears to have a latitudinal contrast with the impact of LAI. While the impact of LAI on the coupling shows strong positive IF values over the NH in the low-end SSPs (figures 4(a3)-(e3)), the impact of SM on the coupling shows positive IF values over the SH in the first period. The historical results show that SM tends to reduce the strength of the meanT-CO 2 causality over the NH and amplifies it over the SH. In SSP1-2.6, SM changes reduce the increases in CO 2 variability due to SM influences on meanT changes (meanT-CO 2 |SM) mainly in the NH in the first period and extend into the SH by the end of the second period of the 21st century (figures 4(b3) and 4(b4)). This could probably be because SM serves as a cooling effect on nearsurface air temperature, which changes its impact on CO 2 variabilities by increasing photosynthesis. Positive IFs increase in the NH during the first period with increasing emission scenarios (figures 4(b3)-(e3)). In the high-end scenarios, the negative IF values over the NH change into positive IF values in the second period of the century. Green et al (2019) suggested that the carbon sinks due to SM might shift toward source functions from mid-century leading to CO 2 growth in the positive IF regions of figures 4(b3)-(e3) possibly owing to the impact on photosynthesis. We find a contrast in these patterns in the SH due to seasonally varying hemispheric conditions which characterize SM (Miralles et al 2012) and vegetation (Wu et al 2015) interaction with the atmosphere.

Cross-comparison of the influences of LAI and SM on the couplings
To highlight what the influences of LAI and SM on the separated couplings look like, their causalities on the CO 2 and meanT couplings are binned over historical LAI and SM climatologies respectively, after which we plot them against each other for the historical (black), SSP1-2.6 (green), and SSP5-8.5 (red) for both projection periods (figure 5). Figures 5(a) and (b) represent the impact on the CO 2 -meanT coupling, and figures 5(c) and (d) represent the impact on the meanT-CO 2 coupling. The results demonstrate a positive linear relationship between them, which intensified in the second century period. SSP1-2.6 and historical results have very similar causal strengths in the first period. However, both factors are observed in the second period to reduce temperature variability amplification due to changes in CO 2 , albeit stronger in SSP1-2.6, shown by the increases in the negative IF values ( figure 5(b)). In the high-end projection, SSP5-8.5, the CO 2 -meanT causality mainly increases. Thus, in the SSP1-2.6, LAI and SM reduce the coupling strength possibly by being sinks to CO 2 , which reduces with increasing emissions (figures 5(a) and (b)). Overall, the results in figures 5(a) and (b) suggest that both factors function similarly to constrain the causality of CO 2 on meanT, which was found to be strongest in SSP1-2.6 and make the CO 2 -meanT in SSP5-8.5 more uncertain with a positive IF rate.
In the meanT-CO 2 causality, we find an inverse relationship with how SM and LAI impact the highand low-end scenarios of the meanT-CO 2 causality, as seen in figures 5(c) and (d). As discussed above and seen in figures 5(a) and (b), the impacts in figures 5(c) and (d) also increase from the first to the second period relative to the historical records. Figure 5 (d) shows that these two factors play significant roles in modulating the coupling at the end of the century, whereby SM intensifies the meanT-CO 2 causality in SSP5-8.5 but LAI reduces it in SSP1-2.6. Contrary to the impact on the CO 2 -meanT coupling, figures 5(c) and (d) show that when LAI serves as a sink, SM serves as a source in the meanT-CO 2 of SSP5-8.5. The opposite is seen for SSP1-2.6.

Conclusion
In this study, we investigated the causal structure between near-surface air temperature and CO 2 using localized spatial scales (pixel-wise) over short time scales (seasonal mean variations of monthly anomalies where their climatologies are removed) using an IF causality approach. Additionally, these localized analyses provided indications of sink-source patterns in the interaction. The projection period (2016-2100) was divided into two: the first period 2016-2060 and the second period 2061-2100 to better understand the evolution of the causal structures over time.
The results indicate that, as identified in previous studies, a mutual causality exists within this interaction, although the spatial and temporal characteristics of the separated causalities are unique. The projection results for the different pathways, representing scenarios of radiative forcing, show that the global CO 2 -meanT causal relation gets stronger with increasing emission scenarios, with the strongest found in SSP5-8.5. Additionally, the areal extent of significant IF from global CO 2 to meanT also increases from just the tropical regions into the high latitudinal regions. Contrasting latitudinal patterns are observed globally in the other direction of the interaction, meanT-CO 2 . These patterns were found to have both spatial and temporal changes for the low-end scenarios (SSP1-2.6, SSP2-4.5) and the high-end scenarios (SSP3-7.0, SSP5-8.5). CO 2 sinks mainly dominated SSP1-2.6 and SSP2-4.5 (negative IF), especially in their second periods where regions that were initially strong positive IFs (CO 2 sources) had become weakly positive or negative IF regions.
On the other hand, SSP3-7.0 and SSP5-8.5 patterns were predominantly strong positive IF values, especially over the NH. To an appreciable degree, the spatial pattern over land agrees with the results obtained in the study by Levy et al (2004), explained as the impact of land use and land cover change on CO 2 variability. This attribution to land cover change was also noted by Zhu et al (2016). Additionally, some studies have reported the possible role of vegetation (here indicated with LAI) and SM on this interaction because of their sink-source functions on CO 2 as well as their impacts on near-surface air temperatures (Zhu et al 2016, Green et al 2019. Thus, we also assessed the roles of these two parameters on the projection outcomes. A detailed study on the impacts of SM and LAI on the interaction revealed that they impact the separated causalities (CO 2 -meanT and meanT-CO 2 ) quite differently. While both function together to constrain the CO 2 -meanT coupling strengths in all the emission scenarios, they inversely affect the meanT-CO 2 coupling strength by the end of the century. This is due to their changing roles as sources and sinks for CO 2 and their modulation of meanT anomalies to either intensify or cool near-surface temperatures reported in previous studies (Zscheischler et al 2015, Schwingshackl et al 2017, Hagan et al 2019. Thus, both SM and LAI can help us understand the resilience of global ecosystems to changes in the CO 2 and meanT interaction within a simplified framework as in this study. However, further studies, like sensitivity experiments, are required to understand the mechanisms leading to these results since these socalled causal tools only assess effects and not mechanisms. Furthermore, it may be necessary to assess how CO 2 emissions during the COVID-19 pandemic would impact these results in future studies as already preliminarily identified by Liu et al (2022). Finally, we note that although Liang (2018) showed that the causality formalism could be applied to a highly nonlinear case to obtain reasonable causal inferences, the results may not be as precise as expected in some cases because of the linear assumptions invoked in the formulation of the approach (Liang 2021).

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// esgf-node.llnl.gov/search/cmip6.