The Many Shades of the Vegetation–Climate Causality: A Multimodel Causal Appreciation

Abstract

The causal relationship between vegetation and temperature serves as a driving factor for global warming in the climate system. However, causal relationships are typically characterized by complex facets, particularly within natural systems, necessitating the ongoing development of robust approaches capable of addressing the challenges inherent in causality analysis. Various causality approaches offer distinct perspectives on understanding causal structures, even when experiments are meticulously designed with a specific target. Here, we use the complex vegetation–climate interaction to demonstrate some of the many facets of causality analysis by applying three different causality frameworks including (i) the kernel Granger causality (KGC), a nonlinear extension of the Granger causality (GC), to understand the nonlinearity in the vegetation–climate causal relationship; (ii) the Peter and Clark momentary conditional independence (PCMCI), which combines the Peter and Clark (PC) algorithm with the momentary conditional independence (MCI) approach to distinguish the feedback and coupling signs in vegetation–climate interaction; and (iii) the Liang–Kleeman information flow (L-K IF), a rigorously formulated causality formalism based on the Liang–Kleeman information flow theory, to reveal the causal influence of vegetation on the evolution of temperature variability. The results attempt to capture a fuller understanding of the causal interaction of leaf area index (LAI) on air temperature (T) during 1981–2018, revealing the characteristics and differences in distinct climatic tipping point regions, particularly in terms of nonlinearity, feedback signals, and variability sources. This study demonstrates that realizing a more holistic causal structure of complex problems like the vegetation–climate interaction benefits from the combined use of multiple models that shed light on different aspects of its causal structure, thus revealing novel insights that are missing when we rely on one single approach. This prompts the need to move toward a multimodel causality analysis that could reduce biases and limitations in causal interpretations.

1. Introduction

Vegetation serves as a vital life support system for human societies, interconnecting soil, atmosphere, and water; however, its interactions within the climate system are very complex [1]. On one hand, global warming profoundly impacts vegetation ecosystems through temporal and spatial changes in climate elements such as water and energy, alongside the ${\mathrm{C}\mathrm{O}}_2$ fertilization effect [2], which contains interactions and potential processes among various influencing factors, contributing to the complexity of the mechanisms by which climate affects vegetation. On the other hand, vegetation ecosystems modulate the climate system through biophysical and chemical processes such as growth, greening, senescence, transpiration, and photosynthesis, thereby amplifying or mitigating signals of climate [3], where the response of vegetation to the climate system is relatively weak and difficult to capture. The interaction between vegetation and climate is thus marked by a complex, bidirectional feedback, with a notable emphasis on the influence of vegetation on temperature dynamics [4]. The predominant effect of climate change on vegetation ecosystems is well documented [5], whereas the reciprocal influence of vegetation on climate change remains underexplored and presents analytical challenges [6]. Methodologically, correlation and regression analysis methods are commonly employed within ecology and geographical sciences to assess the impacts of climate change on vegetation [7,8,9,10], but these methods often falter in clarifying the bidirectional “driving-response” relationship between vegetation and climate, leading to debates and controversies [11]. Climate research, on the other hand, leans toward numerical simulations to decipher vegetation’s feedback to the climate system [12,13], which depend heavily on model parameterization schemes, creating variability in analytical outcomes and obstructing consensus [14].

In light of the aforementioned challenges, causality inferences can either quantitatively and/or qualitatively through interventions unravel cause–effect relationships between various subsystems of the climate [15,16,17]. However, a faithful causal inference is hindered by the complexity of cause–effect relationships, which are characterized by linear and nonlinear relationships, different feedback forms, time- and frequency-dependent components, and sufficient asymmetry needed to untangle the causes from the effects [18]. A good overview of these challenges is presented in [19]. Different causal frameworks have been developed over the years to tackle the different challenges mentioned above and beyond, and while most of them have seen successful applications, each approach is often blind to its limitations. This may indicate the need to move toward a multimodel approach where different approaches with different strengths and limitations are used in tandem to unravel a more comprehensive causal structure.

A classical qualitative causality test approach is the Granger causality (GC) [20], which is a statistical hypothesis test for determining whether one time series is useful in forecasting another. By definition, a variable x (Granger) causes another variable y if the knowledge of y improves the autoregressive forecast of x. Since its development, GC has found a wide application in the fields of finance [21,22,23,24], brain science [25,26,27], and climate science [28,29,30]. For instance, Friston et al. [31] pioneered the detection and estimation of directed connections in neural networks through GC. Kovács et al. [32] employed GC to reveal the underlying causal relationships between environmental drivers and global vegetation attributes, highlighting soil moisture and the availability of accumulated precipitation as critical drivers of vegetation dynamics. Attanasio et al. [33] utilized GC to investigate the relationship between anthropogenic ${\mathrm{C}\mathrm{O}}_2$ emissions and global temperature, uncovering statistical evidence that the increase in ${\mathrm{C}\mathrm{O}}_2$ content leads to a rise in global temperature. Dhamala et al. [34] proposed a time-frequency-based GC to explore multiscale relationship between various variables. However, the assumption of linearity in these earlier forms restricts the scope of Granger causality applications to linear cases. Yet, complex nonlinear feedback relationships often exist in the climate datasets including the vegetation–temperature interaction, prompting the need to develop new methods for inferring nonlinear causal relationships. In that endeavor, Marinazzo et al. [35] applied the theory of reproducing kernel Hilbert spaces to Granger causality analysis to generalize a nonlinear form of GC. This approach was applied to understand the linear and nonlinear characteristics between different variables. Bueso et al. [36] extended kernel Granger causality (KGC) to include timescale applications and applied it to the interactions between El Niño–Southern Oscillation and soil moisture (ENSO-SM).

Land–atmosphere interactions are also characterized by positive and negative feedback signs, which are important for trend analysis in climate signs [37,38]. In general, statistical correlations are often used to represent positive and negative feedback [8]. However, although correlation is a necessary prerequisite for causality, it does not imply causation [39]. Thus, a causal inference model based on a Bayesian network was proposed, by first establishing a Markov graph with correlation to determine the links between network nodes and then through a series of intervention and statistical tests, extract the causal structure from full-time graph, realizing the leap from correlation to causal inference [16]. On this basis, the Peter and Clark momentary conditional independence (PCMCI) method combines the Peter and Clark (PC) algorithm [40] with the momentary conditional independence (MCI) approach to obtain causal graphs [19,41,42,43]. PCMCI has been widely used in many geological and climatological studies [44,45]. For example, Krich et al. [46] utilized PCMCI methods to estimate causal networks in biosphere–atmosphere interactions. Capua et al. [47] used PCMCI to study tropical and mid-latitude teleconnections interacting with Indian summer monsoon rainfall. Qu et al. [48] also employed PCMCI methods to detect the lagged and contingent relationships between wildfire burn area and drought patterns alongside vegetation conditions.

Beyond the above methods, information theoretic (IT) approaches also provide ways to assess causality between time series quantitatively [49]. Schreiber [50] derived an IT measure, the transfer entropy (TE), an asymetric conditional variation of the entropy–based statistics known as mutual information, which has been applied to distinguish drivers from responses within interactions [51,52]. Additionally, IT methods based on entropy, in the Shannon sense [53], provide a way to assess causal histories in systems [54]. Ruddell used mutual information the transfer entropy to assess the cause–effect relationships within a network of ecohydrological variables.

Liang et al. [49] proposed a rigorously formulated causality formalism, from first principles, based on the Liang–Kleeman information flow theory (L-K IF) where causality is quantitatively realized. This has been verified in well-known theoretical models such as the Rossler and Lorenz systems [55], as well its maiden real-world application to reveal the causal relationship between El Nino phenomenon and the Indian Ocean dipole (IOD). Since then, the L-K IF causality has been applied to several areas of research interest in climate science [56], finance [57], and psychology [58], among others, to unravel cause–effect relationships. For instance, Stips et al. [59] used this method to study the causal relationship between greenhouse gas emissions and global warming and found that the increase in the ${\mathrm{C}\mathrm{O}}_2$ concentration in the past 150 years led to the rise of global surface temperature. Hagan et al. [18] applied Kalman filtering and wavelet analysis techniques to extend the Liang–Kleeman information flow method to the time-frequency domain and analyzed the bidirectional mutual feed mechanism between soil moisture and air temperature. Tao et al. [60] used the information flow method to quantitatively assess the effects of climate warming and Atlantic and Pacific decadal oscillations on global precipitation and their regional differences. Docquier et al. [61] applied the Liang–Kleeman information flow method to identify the potential causal drivers of the Arctic sea ice. They found that recent and future changes in Arctic sea ice are primarily driven by air and sea surface temperatures, as well as ocean heat transport, which in turn significantly influence temperature and ocean heat transport dynamics. More recently, Zhou et al. [62] proposed a time-dependent form of the formalism to the multivariate form of the LK causality and applied it to reveal novel findings of time-varying causal structures between soil moisture, vapor pressure deficit, and gross primary productivity.

Understanding natural systems, particularly the climate, is challenging due to the complexity of the systems and the limitations in causal inquiry methods and observations [1]. In scientific research, identifying reliable and relatively certain patterns often relies on extracting and summarizing objective facts [19]. Furthermore, historical observational data provide authentic information and evidence, while causality analysis methods offer rigorous and reliable data analysis and exploration tools [19]. Like many complex natural interactions, the holistic causal structure of the vegetation–climate interaction may not be sufficiently realized with only one approach, since each method comes with its unique strengths and limitations as shown in Figure 1.

In conducting our exploration of the causal relationships between vegetation and climate, we recognize that different causal analysis methods illuminate these causal relationships from varied aspects [63], as shown in Figure 1. This is because causality analysis of natural systems is a very complex problem that has many facets to it. As a result, different causality methods may only provide skewed perspectives of the full causal structure such as the nature of the temporal delays [64,65], degrees of nonlinearity [29], and feedback signs [4,46], among others. Consequently, relying on a single approach limits interpretation of the cause–effect relationship to a specific framework, making it blind to other potentially useful interpretations. This limitation stems from the fact that each method offers a unique lens, focusing on specific aspects of the causal relationship, and thus provides a partial view that falls short of encapsulating the entire complexity of the vegetation–temperature coupling (Figure 1). This study provides a blue–print on how to capture a fuller interpretation of the vegetation–climate coupling by harnessing a multimethod analysis approach.

Nonlinearity is a critical characteristic of the climate system. Therefore, we first conducted a causal analysis of the vegetation–climate system coupling based on the degree of nonlinearity as a function of location. Using the KGC, we aimed to understand global regions of nonlinearity in their causal relationship where knowledge of vegetation improves air temperature prediction. Next, the PCMCI was used to investigate the coupling signs, alongside providing potential temporal delays between the driver and response. Finally, Liang–Kleeman information flow (L-K IF) was applied to highlight regions of uncertainty and predictability in air temperature due to information flow from vegetation changes.

The remaining sections are structured as follows: Section 2 introduces the data and methods applied in this study. Section 3 presents the experimental results. Section 4 discusses and analyzes the experimental findings. Section 5 presents the conclusions.

2. Data and Methods

2.1. Data

2.1.1. CRU TS

For the temperature variable, datasets of global mean air temperature (T) monthly observations at meteorological stations across the world’s land areas, acquired from the Climatic Research Unit (CRU), version TS v4.2.0 [66], were used. CRU TS v4.2.0 Mean Temperature is a high-resolution monthly grid dataset interpolated from 5583 observation stations. It covers the average temperature station data over the last 35 years and was interpolated into 0.5° × 0.5° (longitude × latitude) covering the global land surface. We took the monthly mean temperature from 1981 to 2018 as the temperature variable.

2.1.2. GLASS LAI

For vegetation variables, we utilized LAI data from the Global Land Surface Satellite (GLASS) dataset [67,68]. GLASS LAI employs a generalized regression neural network to train LAI data from satellite observations and AVHRR band data from high-resolution radiometers. The trained model was then applied to estimate the LAI on a global scale. The product closely aligns with mainstream leaf area products, has a spatial resolution of 0.05 degrees and a temporal resolution of 8 days, and spanned the period from 1981 to 2018. In this study, we resampled the spatial resolution of these data to 0.5° × 0.5° (longitude × latitude) to ensure consistency in spatiotemporal resolution with CRU data.

2.1.3. GLEAM ET

Evapotranspiration is very important in land-atmosphere interactions [69]. GLEAM (Global Land Evaporation Amsterdam Model) is a suit of specially designed algorithms aimed at accurately estimating global land evaporation by assimilating soil moisture and vegetation optical depth observations [70,71]. The model divides land evaporation into several key components including transpiration, canopy interception evaporation, bare soil evaporation, ice and snow sublimation and open water evaporation. Here, we employed the transpiration from GLEAM V3.6a to represent evapotranspiration from 1981 to 2018, with a monthly temporal resolution, resampled from a spatial resolution of 0.25° × 0.25° to 0.5° × 0.5° (longitude × latitude), consistent with the vegetation and temperature data specifications.

2.2. Methods

We utilized a series of causality-based analysis methods to construct the global vegetation–temperature causal structure as comprehensively as possible, including KGC, PCMCI, and L-K IF. The main features and differences of these methods are presented in

Table 1. Description of the causal methods used in this paper.

Methods	Kernel Granger Causality	Peter and Clark Momentary Conditional Independence	Liang–Kleeman Information Flow
Abbreviation	KGC	PCMCI	L-K IF
Type of method	Qualitative causality	Qualitative causality	Quantitative causality
Theoretical basis	Granger Causality, spectral representation, kernel function	Conditional independence test, structure causality graph	Liang–Kleeman information flow
Use of time delays	Not by default	Always	Not by default
Use of iterative conditioning	No	Yes	No
Sign meaning 1	No negative value	Positive value: Increases in drivers result in an increase in the target. Negative value: Increases in drivers result in an increase in the target.	Positive value: The driver functions to increase the variability of the target, thereby making it more uncertain. Negative value: The driver functions to reduce variability in the target.
Key references	Marinazzo et al. (2007) [35]	Runge et al. (2019) [19]	Liang (2014) [49]

¹ The sign meaning property of PCMCI and L-K IF is reproduced from [72], with permission from Nonlinear Processes in Geophysics, 2024.

.

The methodological workflow (Figure 2) comprised three steps and is described in the following subsections.

2.2.1. Kernel Granger Causality (KGC)

In GC, a set of stationary time series $\left\{x\right(t){\}}_{t=1,.,N+m}$ can be commonly modeled as [20]

$$X_n=\sum _{i=1}^m{A}_i^{\prime }{X}_{n-i}+\sum _{i=1}^m{B}_i{Y}_{n-i}+E_n^{\prime }$$

(1)

where $Y$ represents another time series that can be similarly defined as $X$. The matrices $A_i^{\prime }$ and $B_i$ represent the autoregression coefficients, while $E_n^{\prime }$ denotes white noise [73]. The parameter m corresponds to the order of the autoregressive model, typically determined through selection criteria such as the Bayesian information criterion (BIC) [74]. A standard least squares optimization method was utilized for estimating model coefficients. Instead of measuring the strength of the causal interaction based on the basic concept of GC, kernel Granger causality (KGC) employs kernel function to project the original linear indivisible time series to reproducing Hilbert spaces ($H$). For each $\alpha \in \{1,\dots ,m\}$, the samples of the $\alpha$th component of $X$ form a vector ${\mathbf{u}}_{\alpha }\in {\mathfrak{R}}^N$. Calling $\mathbf{X}$ the $m\times N$ matrix having vectors ${\mathbf{u}}_{\alpha }$ as rows, $H$ coincides with the range of the $N\times N$ matrix $\mathbf{K}={\mathbf{X}}^{\mathrm{T}}\mathbf{X}$, redefining the Granger causality index as

$${\delta }_{GC}={\displaystyle \frac{{∥P^{\perp }u∥}^2}{1-{\tilde{\mathbf{x}}}^{\mathrm{T}}\tilde{\mathbf{x}}}}$$

(2)

where $\tilde{\mathbf{x}}$ is the projection of an organized vector $\theta ={\left(x_{1+m},\dots ,x_{N+m}\right)}^{\mathrm{T}}$ on Hilbert space $H\subseteq {\mathfrak{R}}^N.$ $\mathbf{P}$ represents the projector on the space $H$. Observing $H^{\perp }$ as the range of the matrix $\tilde{\mathbf{K}}={\mathbf{K}}^{\prime }-\mathbf{P}{\mathbf{K}}^{\prime }-{\mathbf{K}}^{\prime }\mathbf{P}+{\mathbf{P}\mathbf{K}}^{\prime }\mathbf{P}$, the natural choice of the orthonormal basis in $H^{\perp }$ is the set of the eigenvectors of $\tilde{\mathbf{K}}$. Hence, ${∥{\mathbf{P}}^{\perp }\mathbf{x}∥}^2=\sum _{i=1}^m{c}_i^2$, where $c_i$ is the Pearson’s coefficient of $x$ and eigenvectors $t_i$. With a hand of Bonferroni test, a filtered linear Granger causality index can be drawn as

$${\delta }_F(Y\to X)={\displaystyle \frac{\sum _{i^{\prime }}{c}_{i^{\prime }}^2}{1-{\tilde{\mathbf{x}}}^T\tilde{\mathbf{x}}}}$$

(3)

Through spectral representation $K\left(X,X^{\prime }\right)=\sum _a{\lambda }_a{\varphi }_a\left(X\right){\varphi }_a\left(X^{\prime }\right)$, where ${\varphi }_a$ are eigenfunctions of the kernel $k$, it is observed that $\tilde{\mathbf{x}}$, in the feature space spanned by $\sqrt{{\lambda }_a}{\mathrm{\Psi }}_a$, coincides with the nonlinear regression of $x$ in the original variables. Consider the space H spanned by zero-mean vectors $u_{\phi }$, where φ represents eigenvectors of the Gram matrix [70], with elements $K_{ij}=k\left(X_i,X_j\right)$. The Gram matrix $K\prime$ is evaluated using both $X$ and $Y$ to predict $\alpha$, where elements of $K$ are defined as $K_{ij}^{\prime }=k\left(Z_i,Z_j\right)$. The regression values now form a vector that is equal to the projection of α, representing the range of $K\prime$. In subsequent analysis, we consider the inhomogeneous polynomial (IP) of integer order P as our choice to meet our task. The IP kernel [23,24] of integer order $p$ is defined as $K_p\left(X,X^{\prime }\right)={\left(1+X^{\mathrm{T}}{X}^{\prime }\right)}^P$. Along similar lines as described for the linear case, we constructed the kernel Granger causality, taking into account only eigenvectors that passed the Bonferroni test:

$${\delta }_F^K=\sum _{i^{\prime }}{c}_{i^{\prime }}^2$$

(4)

The IP kernel of integer order $p$ is expressed as $X\prime =X \otimes X^T$, and $X^T$ is orthonormal to X (orthogonal unit vectors to $X \otimes X^T$) up to order $p$ and $p$ = 1; consider only the linear regression, where $p \ge 2$ suggests that the information transfer mechanism is nonlinear. The corresponding KGC formula can be expressed as

$$KG{C}_{x\to y}=\sum _t{\left(P^{\perp }y\left(t\right)\right)}^2/\sum _t(y(t)-\overline{y}(\mathrm{t}){)}^2$$

(5)

2.2.2. Peter and Clark Momentary Conditional Independence (PCMCI)

PCMCI combines the Peter and Clark (PC) [40] algorithm momentary conditional independence (MCI) approach. Consider an underlying time-dependent system ${\mathrm{X}}_t=\left(X_t^1,\dots ,X_t^N\right)$ with

$$X_t^j=f_j\left(\mathcal{P}\left(X_t^j\right),{\eta }_t^j\right)$$

(6)

where $f_i$ signifies a potential nonlinear functional dependency, and ${\eta }_t^j$ denotes the mutually independent dynamic noise components; the nodes within a time series graph symbolize the variables $X_t^j$ at various time lags. Here, $\mathcal{P}\left(X_t^j\right)\subset {\mathrm{X}}_t^{-}=\left({\mathrm{X}}_{t-1},{\mathrm{X}}_{t-2},\dots \right)$ is defined to represent the causal progenitors, or “parents”, of the variable $X_t^j$, identified from the historical data across all N variables.

For each variable $X_t^j$, the PC algorithm begins with the establishment of initial parent candidates $\mathcal{P}\left(X_t^j\right)\subset {\mathrm{X}}_t^{-}=\left({\mathrm{X}}_{t-1},{\mathrm{X}}_{t-2},\dots \right)$. Subsequently, the PC algorithm undertakes nonconditional independence evaluations, excising $X_{t-\tau }^i$ if the null hypothesis $X_{t-\tau }^i⫫X_t^j$ stands unrefuted at the significance threshold ${\alpha }_{PC}$. By means of an a-value assessment, causal links in PC algorithms can be defined as

$$p\left(X_{t-\tau }^i\to X_t^j\right)=\underset{\left\{𝒮\right\}}{max}p\left(X_{t-\tau }^i⫫X_t^j\mid 𝒮\right)$$

(7)

The aggregated p-value of a causal link is determined as the maximum of all values resulting from conditional independence tests conducted on various condition sets $𝒮$, as illustrated in Equation (4). Then, the momentary conditional independence (MCI) test is used to test whether $X_{t-\tau }^i\to X_t^j$ with

$$X_{t-\tau }^i⫫/X_t^j\mid \widehat{\mathcal{P}}\left(X_t^j\right)\setminus \left\{X_{t-\tau }^i\right\},\widehat{\mathcal{P}}\left(X_{t-\tau }^i\right)$$

(8)

There are various MCI methods provided by PCMCI, including partial correlation (ParCorr), a nonlinear two-step conditional independence test, and two fully nonparametric tests: (i) a test based on Gaussian process regression and a distance correlation (GPDC) for the detection of additive nonlinear causality and (ii) another test based on conditional mutual information (CMI) for the detection of multiplicative nonlinear causality.

To study the positive and negative feedback relationship between LAI and T, we specifically focused on the linear independence test known as ParCorr.

The ParCorr conditional independence test relies on partial correlations and a t-test, assuming the following model:

$$X^i=𝒮{\beta }_{X^i}+{ϵ}_{X^i},X^j=𝒮{\beta }_{X^j}+{ϵ}_{X^j}$$

(9)

where $\beta$ represents coefficients and ϵ denotes Gaussian noise. This model leads to the subsequent residuals

$$r^{X^i}=X^i-𝒮{\widehat{\beta }}_{X^i},r^{X^j}=X^j-𝒮{\widehat{\beta }}_{X^j},$$

(10)

with $X^i$ being the estimated value. ParCorr eliminates the influence of $𝒮$ on $X^i$ and $X^j$ through ordinary least-squares regression, then testing the independence of the residuals by the Pearson correlation with a t-test. The independence test yields a p-value and a test statistic value I, which represents the correlation coefficient in the case of ParCorr.

2.2.3. Liang–Kleeman Information Flow (L-K IF)

The causality in Liang–Kleeman Information Flow is defined as the time rate of the flow of information transferred from variable $X_2$(e.g., leaf area index (LAI)) to $X_1$ (e.g., temperature (T)) when the system guides a state forward, which can be formulated as:

$${\displaystyle \frac{d{X}_1}{dt}}=F\left(x,t\right)={\displaystyle \frac{d{X}_1^{*}}{dt}}+T_{2\to 1}$$

(11)

where $X_1$ and $X_1^{*}$ are $n$-dimensional vectors, $F={\left(F_1,F_2,\dots ,F_n\right)}^T$ is the vector field.$T_{2\to 1}$ is the information flow (IF) rate from $X_2$ to $X_1$ demonstrated by Liang in [75]. The time rate of change in $X_1$ is precisely equivalent to the mathematical expectation of the divergence of the vector field. Specifically, for the case of bivariate analysis in this study ($n=2$), the information flow (IF) rate is:

$$\begin{array}{c}{T}_{2\to 1}={\displaystyle \frac{d{X}_1}{dt}}-{\displaystyle \frac{d{X}_1^{*}}{dt}} \\ =-E\left({\displaystyle \frac{F_1}{{\sigma }_1}}{\displaystyle \frac{\partial {\rho }_1}{\partial x_1}}\right)-E\left({\displaystyle \frac{\partial F_1}{\partial x_1}}\right)\\ =E\left[{\displaystyle \frac{1}{{\sigma }_1}}{\displaystyle \frac{\partial \left(F_1{\rho }_1\right)}{\partial x_1}}\right] \end{array}$$

(12)

where $E$ denotes the mathematical expectation and ${\sigma }_1$ the marginal probability density of $X_1$. It is important to note that the IF rate is asymmetric $\left(T_{2\to 1}\ne T_{1\to 2}\right)$, distinguishing it from correlation.

In the real cases where the underlying system dynamics are unknown, Equation (2) is estimated using maximum likelihood estimation (MLE). Given two time series $X_1$ and $X_2$, the MLE of $T_{2\to 1}=0$, as described by Liang in [49], is:

$${\widehat{T}}_{2\to 1}={\displaystyle \frac{{\mathrm{C}}_{11}{\mathrm{C}}_{12}{\mathrm{C}}_{2,d1}-{\mathrm{C}}_{12}^2{\mathrm{C}}_{1,d1}}{{\mathrm{C}}_{22}{\mathrm{C}\mathrm{o}}_{11}-{\mathrm{C}}_{12}^2}},$$

(13)

where ${\mathrm{C}}_{ij}$ represents the sample covariance between $X_i$ an $X_j$, and ${\mathrm{C}}_{i,dj}$ is the covariance between $X_i$ and a forward differenced series derived from $X_j$, i.e., ${\dot{X}}_{j,n}=\left(X_{j,n+k}-X_{j,n}\right)/(k\Delta t)$, with $k\ge 1$ as an integer and $\Delta t$ as the time step. Although the Euler forward difference method has lower accuracy, it is often indispensable because it inherently assumes that present conditions can affect future outcomes, but not vice versa.

A nonzero $T_{2\to 1}$ indicates that $X_2$ causes $X_1$. Conversely, $X_2$ is not causal if $T_{2\to 1}=0$, aligning with the principle of nil causality, which is a core theorem in the IF theory, has been empirically tested and verified [55]. As shown in equation (13), if ${\mathrm{C}}_{12}$=0, then $T_{2\to 1}=0$, but the reverse is not true. This highlights the concept that, in a linear sense, causation implies correlation, but correlation does not imply causation [49]. This holds even when the total information flow is zero [76], conclusively addressing the long-standing debate on causation versus correlation. The Fisher information matrix was utilized for significance testing. Its inverse provides a covariance matrix, establishing a significance level at 5%. A significant result indicates that $X_2$ exerts a causal influence on $X_1$ [49].

3. Results

3.1. Nonlinearity of the Causal Relationship between Vegetation and Temperature

The global distribution of the KGC results from LAI to T (LAI→T) from 1981 to 2018 in Figure 3 shows the spatial pattern of vegetation–temperature coupling, presenting differences at various degrees of nonlinearity ((a) P = 1, linear, (b) P = 3, nonlinear, (c) P = 5, strong nonlinear; the rest of the results are provided in Figure S1). Generally, the KGC result of LAI→T is stronger in the mid and low latitudes (intertropical regions) and weaker in the high latitudes. However, stronger signals appear in the high latitudes with increasing nonlinearity (p value). For example, relatively strong KGC values also occur in Eurasia and northern North America as shown in Figure 3b,c.

The linear KGC (P = 1) of the LAI→T distribution shown in Figure 3a predominantly highlights high-value regions in intertropical regions which coincide regions of strong land–atmosphere coupling. These areas were characterized by moderate vegetation cover where changes in LAI could impact changes in T.

In Figure 3b,c, the nonlinear KGC from LAI to T presents a more extensive range of regions with strong LAI responses to temperature, encompassing the intertropical regions, temperate grasslands and parts of high-latitude regions. Notable examples include the central grasslands of North America, temperate forests and grasslands in the Eurasian continent, the Sahel; and southern Africa and southeastern Australia in the Southern Hemisphere. Conversely, regions with weak nonlinear responses of LAI to temperature included extreme arid or cold regions, such as the deserts of Central Asia, and polar regions like Antarctica and Greenland. In these areas, the response of LAI to temperature was minimal, reflected in lower nonlinear GC values. Comparing Figure 3a,c into comparison, it is apparent that the regions of strong LAI→T are characterized by nonlinear processes, which reflects the complexity of the coupling.

For the boreal forest in northern Eurasia, linear KGC analysis at P = 1 (Figure 3a) revealed that the causality strength from LAI to T was comparatively weak. This potentially indicates that changes in T are not directly linked to changes in LAI, and is modulated by other mediating factors which may include precipitation, soil moisture and VPD variability. As a result, P = 5 (Figure 3c) shows stronger signals in the region indicating that the coupling is a potentially complex relationship, and hence, a lower predictive power of T anomalies with LAI.

The box plots in Figure 3d generally show increasing KGC values with increasing P values for LAI→T considering both median and mean distributions. Similar results are also obtained for T→LAI (Figure S1f). These results highlight the need to move beyond linear methods to adequately analyse this interaction. While linear methods may provide some useful information for the nature of the global distribution of the coupling, it may miss out on other key regions which limits their use for regional studies. Fundamentally, these KGC results indicates that accurately predicting T in almost all global regions would require considering more climate factors and processes beyond vegetation since they are mostly nonlinearly linked.

3.2. Feedback/Coupling Signs and Timescales between Vegetation and Temperature

In this section, we compute both the Pearson correlation and PCMCI to represent the linear and nonlinear correlation-based forms respectively. Figure 4 presents the linear Pearson correlation results for the interaction (Figure 4a) and PCMCI results (Figure 4b) for LAI→T globally from 1981 to 2018. We note that Figure 4a, may more rightly represent the feedback relationship while Figure 4b would represent the coupling of LAI to T. Furthermore, the analysis in Figure 4b considered evapotranspiration (ET) as a conditional variable. Positive and negative signals are distinguished by warm red and cool blue colors respectively. White regions are statistically insignificant locations. As shown in Figure 4a, positive feedback signals were predominantly observed in the northern high latitudes, such as Siberia and northern Canada, and the mid-latitude regions, including parts of Europe and North America. Conversely, negative feedback was primarily seen in tropical regions like parts of Brazil, the Sahel and India.

In Figure 4a, in the northern high latitudes and mid-latitude regions, the positive influence of LAI on T suggests that increased vegetation cover led to increases in temperatures, which could be attributed to processes such as reduced albedo, where dense vegetation absorbed more sunlight, leading to increased surface temperatures. Additionally, evapotranspiration from vegetation can contribute to local humidity. During the growing season, increased LAI enhanced photosynthetic activity, leading to higher carbon sequestration and surface cooling through evapotranspiration. This cooling effect can further enhance vegetation growth, creating positive feedback. The coupling observed in these regions, as shown in Figure 4b, may indicate reinforcement of this mechanism, where the delay in response highlighted the time taken for vegetation to significantly influence the local climate. These above potentially illustrate that in the northern high latitudes and mid-latitude regions, energy control was the dominant factor, while water control was relatively weaker. Therefore, the strength of the causal relationship reflected in the water cycle was not as strong as that in the energy cycle. We also note that the spatial distribution of the PCMCI results in agree more with that of the KGC results for P = 5, indicating that the PCMCI does capture the nonlinearity within the interaction.

Conversely, in tropical regions such as Brazil and the Sahelian region, the negative influence of LAI on T was observed. In these areas, increased vegetation cover was associated with lower temperatures. This negative relationship could be explained by the cooling effects of evapotranspiration, where the process of water vapor release from plants led to heat absorption from the surrounding air, thus cooling the surface. Furthermore, dense vegetation could enhance cloud formation, which in turn reflected solar radiation and reduced surface temperatures. These results add a layer of helpful interpretation to the KGC results in Figure 3 which only shows more of a binary causal–noncausal result.

To further analyse the coupling in detail, we selected the northern forest, Central Asian monsoon zone, Sahel region, and Amazon tropical rainforest as representative areas for further analysis. The analysis of four distinct regions—boreal forest (60–65° N, 90–95° E), East Asian monsoon region (26–31° N, 110–115° E), Sahel (5–10° N, 30–35° E), and Amazon rainforest (0–10° S, 55–65° W)—revealed significant insights into the relationship between LAI and temperature T. Across these regions, scatter plots, joint probability density functions (PDFs), and causal structure consistently showed a positive correlation between LAI and T, although with varying degrees of strength and temporal dynamics influenced by regional climatic conditions.

Overall, the plots demonstrate that the interaction is nonlinear, although the degree of nonlinearity varies from region to region. These validate the KGC results. In the boreal forest, Figure 5a shows a robust correlation of 0.88 between LAI and T at a lag of 1 month for LAI→T. Both Figure 5a,b suggest that LAI→T in the region is not a direct or strictly linear relationship. The causal structure (Figure 5c) emphasized a bidirectional causality portrayed as a positive feedback between vegetation and climate in the region.

Similarly, in the East Asian monsoon region, the scatter plots (Figure 5d,e) show a clear positive trend, indicating an even more direct LAI→T relationship at a 1 month lag. Figure 5e shows that this more direct relationship (correlation = 0.96) is likely at all vegetation and temperature scenarios in the region. The causal structure (Figure 5f) here also shows a positive vegetation–climate feedback. In the Sahel, the scatter plots from Figure 5g,h the lagged correlation between LAI and T is −0.69, and the joint PDFs (Figure 5h) demonstrate a multimodal distribution, suggesting complex interactions indicating that multiple environmental factors influences the interaction. The causal structure (Figure 5i) indicated a positive coupling of from T to LAI, and a negative coupling of LAI on regional temperature patterns.

Over the selected Amazon rainforest region, the scatter plots (Figure 5j,k) also highlighted the complexity of vegetation–climate interactions in this region (correlation = −0.43). The PDFs (Figure 5k) also showed a multimodal distribution between vegetation and climate in the region. Like the Sahel, the causal structure (Figure 5l) here also shows positive T→LAI and negative LAI→T. These plots provide a detailed view of the statistical relationship between vegetation and climate which shows the complex relationship between them.

3.3. Information Flow between Vegetation and Temperature

The Liang–Kleeman information flow here was to measure the information flow between LAI and T. Notably, similar to the PCMCI results, the Liang–Kleeman information flow (L-K IF) exhibited both positive and negative values, although their interpretations differ, as indicated in Table 1. L–K IF entropy values qualitative, measured in nats per unit time. As noted in Table 1, for PCMCI, a positive value means that the increasing or decreasing changes in the cause variable results in increases or decreases in the effect, while a negative values denote the contrary. However, in L-K IF, a positive IF indicates that one variable is a source of uncertainty for the other, such that, changes in the causal variable results in increased amplitude or variability of the other variable. Negative IF, on the other hand, suggests that changes in the causal variable reduces the amplitude of the other variable, making the causal variable a source of equilibrium and consequently, predictability. Figure 6 shows the L-K IF from LAI to T from 1981 to 2018 across the globe. As noted above, that the signs of L-K IF, based on entropy, should not be interpreted exactly as that of PCMCI, which is fundamentally based on correlation analysis. Positive values in Figure 6 suggest that vegetation was a source of uncertainty for temperature. Thus, changes in vegetation could result in anomalous temperature events such as heat waves or significant cooling of the temperature in the region. In other words, changes in vegetation could amplify temperature variability in the positive IF regions [77]. Negative L-K IF values suggested that vegetation functions to reduce the variability of temperature, keeping it within a range of equilibrium and predictability as “surprises” or anomalous events in the temperature variability would be less probable (henceforth referred to as the stabilizing effect). Overall, positive values were found in regions of strong land–atmosphere coupling [78,79]. Furthermore, these are also water–limiting regions where changes in soil moisture could potentially drive atmospheric conditions. On the other hand, the blue regions were found in energy-limiting regions across the globe where net radiation drives atmospheric conditions like temperature and precipitation by its control on evaporation [80,81,82]. In Figure 6, northern Eurasia depicts a stabilizing effect of LAI on T, as evidenced by the prevalence of negative L-K IF from LAI to T. This suggests that changes in vegetation in this region contributed to the stabilization of the temperature regime, a phenomenon potentially modulated by evapotranspiration mechanisms (see Figure S4). In northern North America, the L-K IF maps similarly denoted a modest stabilizing role of vegetation on the thermal environment, as indicated by softer blue shades in the LAI→T visual representation. In the Amazon basin, the LAI→T map presented a notably dynamic interaction, with pronounced red-orange gradients which showed that LAI is a source of uncertainty for regional climate system.

4. Discussion

In this study, we used three different causality frameworks to study complex vegetation–climate interactions, demonstrating some of the many facets of causality analysis. Previous studies were based on correlation and regression analysis methods [7,8,9], which remain inherently uncertain. In fact, Li et al. [3] cautioned that it was difficult to conclude that “vegetation greening in northern Eurasia will lead to temperature increase” based on remote sensing data and regression statistical methods. Thus, this study attempts to look into this problem by utilizing causality analysis methods to capture the essential characteristics potentially eliminating spurious causalities based from observational records.

In our study, we took advantage of KGC, which can derive the degree of nonlinearity within an interaction reflecting the different degrees of the vegetation-temperature coupling complexity as a function of location across the globe. The increasing nonlinearity with higher degrees of the polynomial kernel (P, see Figure 3) indicated that the influence of LAI on T involved multifaceted interactions, possibly including other ecological factors and feedback processes. Additionally, our study identified regions of strong nonlinear causal signals consistent with earlier studies like Schwingshackl et al. [83], who found these areas to be characterized by strong soil moisture–evapotranspiration coupling. Combining these findings, we concluded that the hydrological cycle and processes in transitional zones significantly influenced the causal relationship between vegetation and climate, often in a complex, nonlinear manner. The subtleties of these interactions challenge the notion of a straightforward, predictive relationship between vegetation and climate, pointing instead to a system characterized by complex, emergent properties.

Meanwhile, we focused on elucidating the positive and negative feedback and coupling features of the interaction PCMCI, finding that vegetation’s impact on climate exhibited positive feedback and coupling in mainly in the mid- to high-latitude regions, while showing negative feedback in lower mid-to low latitude tropical regions. These results are consistent with recognized global patterns of the feedback between vegetation and air temperature found by Forzieri et al. in [4]. Furthermore, our results showed different feedback signals of vegetation on climate under energy and moisture controls, which is very important in global warming discussions. These findings are also consistent with Krich et al. [46] who also relied on PCMCI. In this study, the combined use of PCMCI and KGC results, provide a detailed analysis of climate transition zones, offering insights into climate change and associated critical scenarios. In regions with strong nonlinearity in land–atmosphere interactions, PCMCI methods extract information on positive and negative feedbacks, along with corresponding time delays, addressing limitations of linear theories and methods that neglect time delays. This detailed analysis, particularly in critical ecological–climatic transition zones globally, enhances our understanding and provides a basis for studying climate change and associated critical scenarios. The L-K IF approach complemented these findings by capturing the dynamics of information flow within the LAI–T coupling. By examining the information flow, this method provided insight into how changes in LAI might induce variability in T, thus indicating regions where changes in LAI could lead to extreme temperature events across the globe, which we find are predominantly in water–limited regions.

In addition to natural factors, human activities also significantly impacted vegetation–climate interactions. However, we note that our study did not take anthropogenic factors into account for the purposes of simplifying the problem. This will be for future studies. The human factor has a great impact on the vegetation on the globe. Real-world scenarios involve complex indirect effects, such as global warming leading to human or animal migration and changes in water use and subsequently affecting urbanization, forest fragmentation, and thus vegetation and ecology [84]. These processes contain important causal information to more adequately understand the causal structure between vegetation and climate. For example, in Figure 4a,b of the study results, there were a few causal signals in the Amazon region that were contrary to the overall regional situation, such as in the southeast of the Amazon in Figure 4a and the scattered positive and negative signals in Figure 4b. These atypical opposite signals may be influenced by human activities, such as the severe human-induced deforestation in the southeastern Amazon found by Ometto et al. in [85] and the impact of road construction on forest fragmentation mentioned by das Neves et al. in [86]. The indirect effects or subprocesses caused by these human factors may be inferred in this paper, but ought to be done with caution.

5. Conclusions

This study systematically studied the causality between vegetation and climate from 1981 to 2018 from observations through different causal characterstics that stem from the added value of different approaches. From these various perspectives, we elucidated the complex impact of vegetation on temperature across the globe. We advocate for an integrated approach to address the climatic and ecological challenges of our time, attaching importance to vegetation dynamics in the climate system. Integrating multiple causality methods, we observed that while KGC highlighted the complexity and nonlinearity in the LAI–T relationship, PCMCI provided a nuanced understanding of the positive or negative feedback, as well as the coupling with time lag, and L-K IF shed light on the system’s equilibrium. The divergence in the conclusions drawn by each method illuminates the multifaceted nature of the LAI–T interaction. The KGC results suggested that the relationship between observed LAI and T is far more complex than a simple cause–effect, depending on the location over the globe. The PCMCI analysis, indicated that the LAI–T relationship was characterized by positive or negative influences. Finally, the L-K IF results on system equilibrium implied that vegetation could be both an agent of uncertainty or equilibrium depending on the location under consideration.

Overall, this study highlights the need to go beyond the absolute conclusions of a single method since each method has its strengths and limitations in providing a targeted perspective on the nature of the causal structure of any system under investigation. Thus, there is the need to employ a multimodel analysis that could bring us closer to obtaining more holistic views of the causal structure of natural systems.