Estimating causal networks in biosphere–atmosphere interaction with the PCMCI approach

The dynamics of biochemical processes in terrestrial ecosystems are tightly coupled to local meteorological conditions. Understanding these interactions is an essential prerequisite for predicting, e.g. the response of the terrestrial carbon cycle to climate change. However, many empirical studies in this field rely on correlative approaches and only very few studies apply causal discovery methods. Here we explore the potential for a recently proposed causal graph discovery algorithm to reconstruct the causal dependency structure underlying biosphere–atmosphere interactions. Using artificial time series with known dependencies that mimic real-world biosphere–atmosphere interactions we address the influence of non-stationarities, i.e. periodicity and heteroscedasticity, on the estimation of causal networks. We then investigate the interpretability of the method in two case studies. Firstly, we analyse three replicated eddy covariance datasets from a Mediterranean ecosystem. Secondly, we explore global Normalised Difference Vegetation Index time series (GIMMS 3g), along with gridded climate data to study large-scale climatic drivers of vegetation greenness. We compare the retrieved causal graphs to simple cross-correlation-based approaches to test whether causal graphs are considerably more informative. Overall, the results confirm the capacity of the causal discovery method to extract time-lagged linear dependencies under realistic settings. For example, we find a complete decoupling of the net ecosystem exchange from meteorological variability during summer in the Mediterranean ecosystem. However, cautious interpretations are needed, as the violation of the method's assumptions due to non-stationarities increases the likelihood to detect false links. Overall, estimating directed biosphere–atmosphere networks helps unravel complex multidirectional process interactions. Other than classical correlative approaches, our findings are constrained to a few meaningful sets of relations, which can be powerful insights for the evaluation of terrestrial ecosystem models.

1 Introduction

The terrestrial biosphere responds to atmospheric drivers such as radiation intensity, temperature, vapour pressure deficit, and composition of trace gases. On the other hand, the biosphere influences the atmosphere via partitioning incoming net radiation into sensible, latent, and ground heat fluxes, as well as via controlling the exchange of trace gases and volatile organic compounds. Over the past few decades, many of these processes have been identified and their physical, chemical, and biological effects have been investigated (see, e.g. Monson and Baldocchi, 2014; McPherson, 2007, for overviews). However, there are still substantial unknowns regarding the exact causal dependencies among the different processes (Baldocchi et al., 2016; Miralles et al., 2018), which leads to large uncertainties when predicting, e.g. ecosystem responses to drought conditions (von Buttlar et al., 2018; Sippel et al., 2017).

Multiple ecological monitoring systems have been set up to monitor ecosystem dynamics. Networks of eddy covariance towers continuously monitor carbon, water, and energy fluxes in high temporal resolution (Baldocchi, 2014). Satellite remote sensing data complement this picture and can be used in tandem (Papale et al., 2015; Mahecha et al., 2017). They typically only monitor vegetation states at multi-day resolutions and some products offer nearly complete global coverage (Justice et al., 2002; Woodcock et al., 2008). The actual and future satellite missions are leading to rapid development in the field with ever higher spatial, temporal, and spectral measurements (Malenovský et al., 2012; Guanter et al., 2015; Qi and Dubayah, 2016).

The study of biosphere–atmosphere interactions using observations typically relies on correlative approaches, or is based on model data, i.e. requires a priori knowledge. In recent years, a new branch in statistics aiming for causal inference from empirical data has experienced substantial progress. The idea of causal inference had already emerged in the early 20th century (Wright, 1921). Later, Granger (1969) suggested one of the first applicable formalisms; since then, several efforts in ecology and climate science have concentrated on the bivariate form of Granger causality (Elsner, 2006, 2007; Kodra et al., 2011; Attanasio, 2012; Attanasio et al., 2012). From an information theory perspective, transfer entropy (Schreiber, 2000) evolved as a frequently used measure to infer directionality and amount of information flow (Kumar and Ruddell, 2010; Ruddell et al., 2015; Gerken et al., 2018; Yu et al., 2019). For instance, Ruddell and Kumar (2009) used transfer entropy to estimate networks of information flow. These networks constructed for an agricultural site under drought and non-drought conditions showed substantial differences in connectivity, especially between subsystems comprised of variables of land and atmospheric conditions. Those changes in connectivity are attributed to changes in the feedback patterns between the subsystems for drought and normal conditions. The original forms of both Granger causality and transfer entropy are bivariate and converge for the case of vector auto-regressive models. While Granger causality is typically limited to linear relationships, transfer entropy also captures non-linear interactions but requires very large data quantities for the estimation of the probability density function.

Aiming to mitigate some of the limitations of the traditional Granger causality, Detto et al. (2012) used a conditional spectral Granger causality framework that allowed them to disentangle system inherent periodic couplings from external forcing. The disentanglement is enabled via decomposition into the frequency domain using wavelet theory. This method led them to the finding that soil respiration in a pine and hardwood forested ecosystem in winter is not influenced by canopy assimilation but only by temperature, a result that would not be detectable via lagged correlation or bivariate Granger causality. A time frequency representation of Granger causality was presented by Shadaydeh et al. (2019), which allowed them to identify anomalous events in marine and ecological time series. Green et al. (2017) used a similar approach to Detto et al. (2012) to investigate biosphere–atmosphere feedback loops. It was found that they can explain up to 30 % of variance in radiation and precipitation in certain regions. Recently, Papagiannopoulou et al. (2017a) applied a non-linear multivariate conditional Granger causality framework to study climatic drivers of vegetation at the global scale. This approach revealed that water availability dominates plant productivity, as 61 % of the vegetated surface appeared water-limited rather than controlled by radiation or temperature (Papagiannopoulou et al., 2017b). In cases of transfer entropy, Goodwell and Kumar (2017a, b) developed a redundancy measure, which allows them to distinguish unique, synergistic, and redundant information transfer of a bivariate or (potentially) multivariate system to a target variable. This modification enables a stronger multivariate interpretation of process networks constructed with transfer entropy. Changes in connectivity then potentially point to different ecosystem response strategies to disturbances (Goodwell et al., 2018). These examples highlight that unexpected interaction patterns can in principle be identified from data only and may challenge theoretical assumptions. In fact, in the last few years the science of causal inference has developed a strong theoretical foundation and several algorithms have been proposed (Spirtes et al., 2001; Pearl, 2009; Peters et al., 2017; Pearl and Mackenzie, 2018; Runge et al., 2019a). However, only a few studies test the suitability of this latest generation of methods for understanding ecosystem dynamics (see, e.g. Shadaydeh et al., 2018, 2019; Christiansen and Peters, 2020).

Ecological and climate data are often time-ordered. This property can be exploited to construct time series graphs (Ebert-Uphoff and Deng, 2012). Recently, Runge et al. (2019b) introduced an algorithm to estimate such graphs, called PCMCI, a combination of the PC algorithm (named after its inventors Peter and Clark; Spirtes and Glymour, 1991) and the Momentary Conditional Independence (MCI) test. PCMCI has been successfully applied to artificial (Runge et al., 2018) and climatological case studies (Runge et al., 2014; Kretschmer et al., 2016). Hence, this method could be potentially of high relevance for learning the causal dependency structure underlying biosphere–atmosphere interactions.

In this study, we explore the potential of PCMCI for disentangling and quantifying interactions and feedbacks between terrestrial biosphere state and fluxes and meteorological variables. The study is structured as follows. In Sect. 2, we justify and introduce the method from an ecological perspective, and we also describe artificial and real-world datasets explored in this study. The results in Sect. 3 describe the performance of the method on artificial time series data with known dependencies that mimic some basic properties of observed land surface fluxes such as heteroscedasticity. We then report on the exploration of three replicated eddy covariance measurement towers in a Mediterranean ecosystem and explore how the identified interdependencies of carbon and energy fluxes and micro-meteorological observations vary over time. Further, we present the analysis of global satellite data of vegetation greenness to understand the lagged dependency of ecosystems with respect to climatic drivers. Based on these results, Sect. 4 discusses the potentials and limitations of PCMCI for other applications in land–atmosphere studies and give recommendations for further methodological developments.

2 Method and data

2.1 From bivariate to multivariate measures of causality

Monitoring an ecosystem with continuous observations of net ecosystem exchange (NEE), the underlying gross primary production (GPP), and ecosystem respiration (R_eco) together with the relevant drivers, i.e. global radiation (R_g), surface air temperature (T), and soil moisture (SM), allows us to study the dynamics of the carbon cycle in terrestrial ecosystems. To foster its understanding, a fundamental question is how these variables causally depend on each other. This requires the identification of directional dependencies such as the well-known effects of SM → GPP and GPP → R_eco and their differentiation from physically implausible links such as R_eco → GPP. Graphical causal models (Spirtes et al., 2001) provide a framework to represent and identify causal relations based on conditional independence relations in data streams of this kind. In the case of an ecological monitoring site as described here, we can exploit the temporal information of the observations for identifying a time series graph as a type of graphical model (Runge, 2018). Formally this can be stated as follows. The variable $X_t^i$ is part of a multivariate stochastic process X (where i is the variable index, in the example $i\in \mathit{\{}{R}_{\mathrm{g}},T,\mathrm{GPP},R_{\mathrm{eco}},\mathrm{SM}\mathit{\}}$, and t is the time index). A time series graph 𝒢 visualises how the individual variables Xⁱ∈X depend on each other at specific time lags τ, i.e. $X_{t-\mathit{\tau }}^i$ with $\mathit{\tau }\in \mathit{\{}\mathrm{1},\mathrm{\dots },{\mathit{\tau }}_{max}\mathit{\}}$ (see Runge, 2018, for definitions). In the following, we refer to a variable $X_{t-\mathit{\tau }}^i$ that is causally affecting a variable $X_t^j$ as “parent” or “driver” and the latter as “receiver” or “target”. To come to a causal interpretation, it is important to exclude dependencies between two variables that are due to common drivers ($X_{t-{\mathit{\tau }}_{\mathrm{1}}}^i\leftarrow X_{t-{\mathit{\tau }}_{\mathrm{2}}}^s\to X_t^j$) or indirect paths ($X_{t-{\mathit{\tau }}_{\mathrm{2}}}^i\to X_{t-{\mathit{\tau }}_{\mathrm{1}}}^s\to X_t^j$). For instance, when estimating the effects of GPP on R_eco and R_g on R_eco using a bivariate measure, one likely obtains implausible results like links that are too strong or even unexpected links because T, respectively, as the common driver and mediator (indirect path), is not accounted for. To exclude dependencies due to common drivers or indirect paths, conditional independence tests are used, denoted as $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j|\mathcal{S})$, with some conditioning set 𝒮. If any variable (or their combination) in 𝒮 explains the dependence between $X_{t-\mathit{\tau }}^i$ and $X_t^j$, then the CI statistic is 0.

Two prominent methods that aim for directional dependencies are Granger causality and transfer entropy (Granger, 1969; Schreiber, 2000). Granger causality is typically estimated as a vector auto-regressive model and thus captures only linear links. Transfer entropy, based on information theory, also captures non-linear dependencies. It can be shown that for multivariate Gaussians, transfer entropy is equivalent to Granger causality (Barnett et al., 2009). Both can be phrased as testing for conditional independence (Runge et al., 2019a). In their original bivariate form, neither of these two methods accounts for third variables, but both can also be extended to deal with multivariate time series as required here (Runge et al., 2012; Granger, 1969). There are even non-linear and spectral modifications of Granger causality which have been applied to study biosphere–atmosphere interactions (Papagiannopoulou et al., 2017a; Detto et al., 2012; Claessen et al., 2019). However, the estimation of multivariate transfer entropy is challenging due to the “curse of dimensionality” (Runge et al., 2012), and multivariate Granger causality also exhibits low link detection power for larger number of variables (higher dimensions) and limited sample size, as is the case in our application (Runge et al., 2019b). The strong decrease in detection power happens when using the whole past ${\mathbf{X}}_t^{-}=({\mathbf{X}}_{t-\mathrm{1}},{\mathbf{X}}_{t-\mathrm{2}},\mathrm{\dots })$ of $X_t^j$, truncated at a maximum lag τ_max, as a conditioning set 𝒮. The problem is that this set can contain a high number of conditions that are irrelevant. For example, when assessing the effect of R_g at a specific time lag τ on GPP using multivariate Granger causality one would create a vector auto-regressive model comprised of all variables, i.e. R_g, T, SM, GPP, and R_eco, at each available lag, but R_eco, dominated by heterotrophic respiration, is not expected to affect gross primary productivity and could be removed to decrease the dimensionality. However, manually selecting conditions is not desirable when the underlying dependence structure is unknown, which is why ideally the conditioning set is identified automatically.

PCMCI addresses this issue by reducing the set of conditions 𝒮 prior to quantifying the dependence between two variables. The two-step approach utilises a variant of the PC algorithm (Spirtes and Glymour, 1991) and the momentary conditional independence measure (MCI) (Runge et al., 2019b). More detailed descriptions are given in Sect. 2.4 and 2.5, respectively (a full description of PCMCI including proofs and quantitative comparisons with other methods are provided in Runge et al., 2019b). A schematic of the PCMCI approach is given in Scheme A1 in Appendix A. PCMCI belongs to the family of causal graphical models (Spirtes et al., 2001; Pearl, 2009) and follows the assumptions listed in Sect. 2.2. In the limit of infinite time series length, PCMCI converges to the true graph of dependencies, which is why we use the term “causal”. As we deal with finite sample length and partially unfulfilled assumptions, spurious links can still appear (beyond the expected false positive rate) and therefore each detected link has to be interpreted with caution.

2.2 Assumptions

PCMCI assesses the causal structure of a multivariate dataset or process X by estimating its time series graph. To draw causal conclusions from observational data, any causal method must adopt a number of assumptions (Pearl, 2009; Spirtes et al., 2001). For the time series case, here we assume time order, the causal Markov condition, faithfulness, causal sufficiency, causal stationarity, and no contemporaneous causal effects (Runge et al., 2018). PCMCI is applied in combination with the ParCorr linear independence test based on partial correlations (see Sect. 2.3). This application additionally requires stationarity in the mean and variance and linear dependencies. In the following, we briefly discuss these assumptions (further details in Runge, 2018; Runge et al., 2019b).

The time order within the time series allows us to orient directed links which are only pointing forward in time. This accounts for causal information propagating forward in time only, i.e. the cause shall precede the effect. Therefore, a directed causal link $X_{t_{\mathrm{1}}}^i\to X_{t_{\mathrm{2}}}^j$ can only exist between two nodes $X_{t_{\mathrm{1}}}^i,X_{t_{\mathrm{2}}}^j$ if t₁<t₂. When a contemporaneous link is found, i.e. t₁=t₂, it is considered to be undirected. In ecological language it means that in order to claim that R_g is driving GPP, any change in R_g that is affecting GPP must be measured at a time step before the change in GPP occurs. The causal Markov and faithfulness assumptions relate the underlying physical causal mechanisms to statistical relationships manifest in the data. The causal Markov condition states that if two processes are not directly connected by some physical mechanism, then they should be statistically independent conditional on their direct drivers, like R_g and R_eco conditional on T. The faithfulness assumption concerns the other direction: if two processes are statistically independent, then there cannot be a direct physical mechanism. The causal sufficiency assumption implies that every common cause of two or more variables Xⁱ∈X is included in X. If this is not the case, detected links may be indirect or due to an unobserved common driver. However, the absence of a link in the detected graph still implies that no direct link is present (as this only requires the assumption of faithfulness). For example, R_g is expected to influence R_eco via T, the indirect path. However, a link between R_g and R_eco might be detected if T is not included in the analysis. However, a missing link between T and R_eco might indicate conditions inhibiting respiratory processes, i.e. very cold temperatures with frozen surfaces or very dry conditions with dead vegetation coverage. Causal stationarity refers to the existence of links over time. In a deciduous forest, for example, the ecosystem's CO₂ exchange is not causally stationary as the link R_g → GPP is given in summer but not in winter. Formally, a process X with graph 𝒢 is called causally stationary over a time index 𝒯 if and only if for all links $X_{t-\mathit{\tau }}^i\to X_t^j$ in the graph the condition $X_{t-\mathit{\tau }}^i\textcolor{red}{}{X}_t^j|{\mathbf{X}}_t^{-}\backslash \mathit{\{}{X}_{t-\mathit{\tau }}^i\mathit{\}}$ holds for all t∈𝒯.

2.3 Independence test

At the core of PCMCI there are conditional independence tests $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$ to evaluate whether $X_{t-\mathit{\tau }}^i\textcolor{red}{}{X}_t^j|\mathcal{S}$ given a conditioning set 𝒮. Within the PCMCI software package TIGRAMITE (Runge, 2017), several independence tests are implemented. Here, we focus on the linear independence test called ParCorr. The ParCorr conditional independence test is based on partial correlations and a t test. This assumes the following model:

$$\begin{array}{}\text{(1)}& X^i=\mathcal{S}{\mathit{\beta }}_{X^i}+{\mathit{ϵ}}_{X^i},X^j=\mathcal{S}{\mathit{\beta }}_{X^j}+{\mathit{ϵ}}_{X^j},\end{array}$$

with coefficients β and Gaussian noise ϵ. This leads to the following residuals:

$$\begin{array}{}\text{(2)}& r^{X^i}=X^i-\mathcal{S}{\widehat{\mathit{\beta }}}_{X^i},r^{X^j}=X^j-\mathcal{S}{\widehat{\mathit{\beta }}}_{X^j},\end{array}$$

with estimated $\widehat{\mathit{\beta }}$. ParCorr removes the influence of 𝒮 on Xⁱ and X^j via ordinary least-squares regression and tests for independence of the residuals using the Pearson correlation with a t test. The independence test returns a p value and test statistic value I, i.e. the correlation coefficient in case of ParCorr. Thus, to identify the effect of GPP on R_eco that does account for their common driver T∈𝒮, ParCorr will perform a linear regression of T on both GPP and R_eco accounting for time lags. The p value of the residuals' partial correlation test can be used to assess whether the two variables are dependent.

2.4 PC algorithm

To efficiently estimate $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j|\mathcal{S})$ the conditioning set 𝒮 should be as small as possible, which means that it should only contain relevant conditions, which allows us to isolate the unique influence of $X_{t-\mathit{\tau }}^i$ on $X_t^j$. For an estimation of $\mathrm{CI}({R_{\mathrm{g}}}_{t-\mathit{\tau }},{\mathrm{GPP}}_t|\mathcal{S})$, for example, 𝒮 should contain T and SM (at certain lags), as they influence the ability of an ecosystem to perform photosynthesis. Likewise, when estimating $\mathrm{CI}(T_{t-\mathit{\tau }},{\mathrm{GPP}}_t|\mathcal{S})$, 𝒮 should include R_g and SM for the same reasons. A sufficient set of relevant conditions includes the drivers and parents of the variable $X_t^j$. Consequently, the aim of the PC step is to identify an as small as possible superset of the parents of each variable included in the process. The algorithm uses a variant of the PC algorithm (Spirtes et al., 2001); a comprehensive pseudo-code of this procedure is given in the Supplement of Runge et al. (2019b). In the limit of infinite sample size, the relevant conditions indeed converge to the true causal parents, although, practically, an estimate that contains a few irrelevant conditions, like R_eco, is sufficient as well.

The PC step starts by initialising the whole past of a process: $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)={\mathbf{X}}_t^{-}=\mathit{\{}{X}_{t-\mathit{\tau }}^i:i=\mathrm{1},\mathrm{\dots },N,\mathit{\tau }=\mathrm{1},\mathrm{\dots },{\mathit{\tau }}_{max}\mathit{\}}$. Next, by evaluating $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$, conditions $X_{t-\mathit{\tau }}^i$ are removed from $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ that are independent of $X_t^j$ conditionally on a subset $\mathcal{S}\in \stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)\backslash \left\{X_{t-\mathit{\tau }}^i\right\}$. 𝒮 starts as the empty set ∅ and is iteratively increased. For instance, let $X_{t-\mathit{\tau }}^i$ be R_eco (at a specific lag) and $X_t^j$ be GPP. The conditional independence between GPP and R_eco will be estimated first by using no conditions. If GPP and R_eco appear related, one variable will be included in the conditioning set. If the residuals are still dependent, a second variable is included and so on. When T is part of the conditioning set, the residuals of GPP and R_eco might not be dependent anymore and R_eco is removed from the estimated set of parents of GPP. The PC algorithm adopted in PCMCI efficiently selects those conditioning sets to limit the number of tests conducted.

Every conditional independence test is evaluated at a significance threshold α_pc, which is usually set to a liberal value between 0.1 and 0.4. Alternatively, in TIGRAMITE one can let α_pc be unspecified. PCMCI then evaluates the best choice of α_pc $\in \mathit{\{}\mathrm{0.1},\mathrm{0.2},\mathrm{0.3},\mathrm{0.4}\mathit{\}}$ based on the Akaike information criterion (AIC), which is further explained in Runge et al. (2018).

2.5 MCI tests

MCI is the actual causal discovery step that ascribes a p value and strength to each possible link. MCI iterates through all pairs $(X_{t-\mathit{\tau }}^i,X_t^j),i=\mathrm{1},\mathrm{\dots },N,\mathit{\tau }=\mathrm{0},\mathrm{\dots },{\mathit{\tau }}_{max}$, and calculates $\mathrm{CI}(X_{t-\mathit{\tau }}^i,X_t^j,\mathcal{S})$, where 𝒮 consists of the two (super-)sets of parents $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ and $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_{t-\mathit{\tau }}^i\right)$ obtained in the PC step. $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_{t-\mathit{\tau }}^i\right)$ is constructed by shifting the time series of $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^i\right)$ by τ. In cases where $X_{t-\mathit{\tau }}^i\in \stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$, $X_{t-\mathit{\tau }}^i$ has to be removed from $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$. If τ=0, conditional dependence is estimated for contemporaneous nodes $X_t^j$ and $X_t^i$. Due to missing time order, a dependence would be left undirected. Further, as the parents $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^i\right)$ and $\stackrel{\mathrm{̃}}{\mathcal{P}}\left(X_t^j\right)$ used in each conditional dependence test are defined to lie in the past of $X_t^i$ and $X_t^j$, links, both contemporaneous and lagged, can be spurious due to contemporaneous common drivers or contemporaneous indirect paths. The absence of a link, though, means that a physical (contemporaneous) link is unlikely (assuming faithfulness; see Runge et al., 2018). For simplicity, the previously given examples were omitting the time lag. Thus, if R_eco responds instantaneously (considering the sampling temporal resolution) to changes in T but T responds with a time lag to R_g, both variables will likely appear contemporaneously coupled to R_g.

The link strength in the PCMCI framework can be given by the effect size of the conditional independence test statistic measure CI used in combination with MCI. In case of ParCorr, the effect size is given by the partial correlation value, which is between −1 and 1. Assuming a linear Gaussian model the partial correlation value is shown to directly depend on the receiver's and driver's variance as well as the coupling coefficient (Runge et al., 2019a):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\rho }}_{X_{t-\mathit{\tau }}^i\to X_t^j}^{\mathrm{MCI}}{\displaystyle }={\displaystyle \frac{c{\mathit{\sigma }}_{X_{t-\mathit{\tau }}^i}}{\sqrt{{\mathit{\sigma }}_{X_t^j}^{\mathrm{2}}+c^{\mathrm{2}}{\mathit{\sigma }}_{X_{t-\mathit{\tau }}^i}^{\mathrm{2}}}}},\end{array}$$

where ${\mathit{\sigma }}_{X_{t-\mathit{\tau }}^i,X_t^j}$ are the variances of the noise terms driving $X_{t-\mathit{\tau }}^i$ and $X_t^j$, respectively, and c is their coupling coefficient. In practice, non-linear links can often also be well detected with ParCorr in so far as they are linearly approximated. In case the linear part is even “stronger” than the non-linear part, ParCorr might also have a better detection power than a non-linear independence test (Runge et al., 2018).

2.6 Data

2.6.1 Artificial dataset – test model

We tested the algorithm on artificial datasets prior to its application to real-world data. The artificial dataset was created using a test model that takes time series of measured global radiation (ℛ_g) and creates three artificial time series that conceptually represent temperature (𝒯), gross primary production (𝒢𝒫𝒫), and ecosystem respiration (ℛ_eco). Note that this test model is not intended to accurately represent observed land–atmosphere fluxes but only serves to test the procedure. The model incorporates one linear auto-dependence ${\mathcal{T}}_{t-{\mathit{\tau }}_{\mathrm{1}}}\to {\mathcal{T}}_t$, one linear additive cross-dependence ${{\mathcal{R}}_{\mathrm{g}}}_{t-{\mathit{\tau }}_{\mathrm{2}}}\to {\mathcal{T}}_t$, and two non-linear dependencies, multiplicative ${{\mathcal{R}}_{\mathrm{g}}}_{t-{\mathit{\tau }}_{\mathrm{3}}}\cdot {\mathcal{T}}_{t-{\mathit{\tau }}_{\mathrm{4}}}\to \mathcal{G}\mathcal{P}{\mathcal{P}}_t$ and multiplicative exponential $\mathcal{G}\mathcal{P}{\mathcal{P}}_{t-{\mathit{\tau }}_{\mathrm{5}}}\cdot c^{{\mathcal{T}}_{t-{\mathit{\tau }}_{\mathrm{6}}}}\to {{\mathcal{R}}_{\mathrm{eco}}}_t$ (see Fig. 1), according to the following equations:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathcal{R}}_{{\mathrm{g}}_{\mathrm{mo}}}& {\displaystyle }={\mathcal{R}}_{{\mathrm{g}}_{\mathrm{obs}}},\text{(5)}& {\displaystyle }{\mathcal{T}}_{\mathrm{mo}}\left(t\right)& {\displaystyle }=c_{\mathrm{1}}{\mathcal{T}}_{\mathrm{mo}}(t-{\mathit{\tau }}_{\mathrm{1}})+c_{\mathrm{2}}{\mathcal{R}}_{{\mathrm{g}}_{\mathrm{mo}}}(t-{\mathit{\tau }}_{\mathrm{2}})+{\mathit{\xi }}_T,\text{(6)}& {\displaystyle }\mathcal{G}\mathcal{P}{\mathcal{P}}_{\mathrm{mo}}\left(t\right)& {\displaystyle }=c_{\mathrm{3}}{\mathcal{R}}_{{\mathrm{g}}_{\mathrm{mo}}}(t-{\mathit{\tau }}_{\mathrm{3}})\cdot {\mathcal{T}}_{\mathrm{mo}}(t-{\mathit{\tau }}_{\mathrm{4}})+{\mathit{\xi }}_{\mathrm{GPP}},\text{(7)}& {\displaystyle }{\mathcal{R}}_{{\mathrm{eco}}_{\mathrm{mo}}}\left(t\right)& {\displaystyle }=c_{\mathrm{4}}\mathcal{G}\mathcal{P}{\mathcal{P}}_{\mathrm{mo}}(t-{\mathit{\tau }}_{\mathrm{5}})\cdot c_{\mathrm{5}}^{\frac{{\mathcal{T}}_{\mathrm{mo}}(t-{\mathit{\tau }}_{\mathrm{6}})-T_{\mathrm{ref}}}{\mathrm{10}}}+{\mathit{\xi }}_{{\mathcal{R}}_{\mathrm{eco}}}.\end{array}$$

The parameters c₁, c₂, …, c₅ are referred to as coupling coefficients, and the time lags are noted as τ₁, τ₂, …, τ₅. The subscripts mo and obs abbreviate model and observation, respectively. T_ref is set to 15 ^∘C. The term ξ, termed “intrinsic” or “dynamical noise”, here represents values from uncorrelated, normally distributed noise. Having dynamic noise is essential for a method utilising conditional independence tests. It is based on the assumption that a process or state is never fully described by its deterministic part because there are unresolved intrinsic processes, summarised as ξ.

The model was fitted to real observational data (using radiation, temperature, and land–atmosphere fluxes) of daily time resolution, measured by the eddy-covariance method (Baldocchi et al., 1988; Baldocchi, 2003) from FLUXNET, by minimising the sum of squared residuals using the gradient descent implemented in the Optim.jl package (Mogensen and Riseth, 2018). We fitted the model to 72 sites listed in Table A2 given in Appendix A. The value range for the coupling coefficients c₁ to c₄ and c₅ were set to [0.2,1] and [1,2.5], respectively. The lags were limited to integer values in the range [0,25]. The distributions of obtained lags and coupling coefficients are given in Figs. A2 to A5 in Appendix A. The fitting process thus generated 72 sets of parameters, containing coupling coefficients and lag values, which were used for the time series generation.

From each of the 72 sets of parameters we generated four sets of time series, each having a length of 500 years. The time series generation was initiated using two types of data: first, uncorrelated, normally distributed noise, and second, unprocessed radiation data as used during the fitting (the available radiation data were repeated to 500 years). The resulting datasets are called baseline dataset and seasonality dataset, respectively. In both cases, the model was run twice, once with homoscedastic (constant variance of ξ), once with heteroscedastic dynamical noise ξ. To induce heteroscedasticity, ξ was multiplied with a mean daily variance that was extracted for each variable at each FLUXNET site. In Fig. A6, a 5-year time series excerpt from the Hainich site (Knohl et al., 2003b) is shown. A third dataset is generated by anomalisation (subtraction of smoothed seasonal mean) of the seasonality dataset.

2.6.2 Eddy covariance data – Majadas de Tiètar experimental site

Data from three towers located in Majadas de Tiètar, (ES-LMa, ES-LM1, ES-LM2), a Mediterranean Savanna in central Spain, are used (coordinates of the central tower: 39^∘56^′25^′′ N, 5^∘46^′29^′′ W). Measurements include the exchange of CO₂ between the land surface and atmosphere at half-hourly resolution using the eddy covariance method. The three tower footprints received different fertilisation treatments in spring 2015 (El-Madany et al., 2018). We consider data from before the fertilisation from April 2014 to March 2015 of shortwave downward radiation (R_g), air temperature (T), net ecosystem exchange (NEE), vapour pressure deficit (VPD), sensible heat (H) and latent heat (LE). The average temperature within this period was 17.3^∘ C, with a total precipitation of 765 mm. Most precipitation fell between October and April.

We expect the causal imprints in the data to vary between seasons and during the course of the day. To satisfy causal stationarity, we estimate networks separately for each month and consider only samples for which the potential radiation was above four-fifths of the potential daily maximum, which corresponds to midday samples. We used a mask type that limits only the receiver variable to the respective month and day time values (see Table A1 for PCMCI parameter settings). This setting causes time series lengths ranging from 239 datapoints in December to 372 datapoints in July. Minimal and maximal lags were set to 0 and 8, respectively. We left the data unprocessed, i.e. we did not subtract a seasonal mean for anomalisation. Constraining the samples to separate month and midday values reduces the effect of seasonality as a common driver that would lead to spurious links. Furthermore, to correct for multiple testing we applied the Benjamini–Hochberg false discovery rate correction (Benjamini and Hochberg, 1995). Thereby, the p values for the whole graph obtained from the MCI step are adjusted to control the number of false discoveries (Runge et al., 2018). We chose a two-sided significance level of 0.01.

2.6.3 Gridded global dataset

The second observational case study was performed on a global dataset. We used data with 0.5^∘ spatial and monthly temporal resolution from 1982 to 2008. The dataset is composed of three climatic variables, global radiation (R_g), temperature (T) and precipitation (P), and one vegetation state index, the Normalised Difference Vegetation Index (NDVI). Both temperature and precipitation datasets were taken from the Climate Research Unit (CRU), version TS3.10 (Harris et al., 2014). The radiation data stems from the Climate Research Unit and National Centers for Environmental Prediction dataset (CRUNCEP, Viovy, 2016). The used NDVI data stems from the Global Inventory Monitoring and Modeling Systems (GIMMS) in version 3g_v1 (Pinzon and Tucker, 2014).

To examine the influence of radiation, temperature, and precipitation on NDVI by means of PCMCI, we used the following settings. We compute the anomalies by subtracting a smoothed seasonal mean. A maximal time lag of 3 months was chosen based on the largest lag with significant partial correlation among all pairs of variables, partialling out only the autocorrelation of each variable. The receiver variable was limited to the growing season defined by T>0 and NDVI>0.2, which allows good comparison to Wu et al. (2015). The significance level (α_pc) in the condition selection phase (see Sect. 2.4) was chosen based on the AIC selection criterion. A concise list of PCMCI parameters that were altered from default settings is given in Table A1.

3 Results

3.1 Test model

As an example, in Fig. 1 we show PCMCI and lagged correlation networks in the form of process graphs. It is clearly visible that many more spurious links pass the significance threshold of 0.01 using lagged correlation as compared to using PCMCI. Those spurious links can complicate the analysis or lead to false conclusions and misleading hypotheses. We examined four cases of different time series lengths: 91 [$\mathrm{183}\le \mathrm{doy}\le \mathrm{274}$] and 120 [$\mathrm{153}\le \mathrm{doy}\le \mathrm{274}$] d and 1 and 5 years for daily data (doy: day of the year). For each time series length and each parameter set, the causal network structure was estimated for 100 realisations of the model (each based on a realisation of intrinsic noise), which allowed the estimation of false positive (FPR) and true positive (TPR) detection rates. The detection rates are calculated for each tower, FPR in general and TPR link-wise. The TPR for each link is its sum of detections among 100 realisations divided by 100. The FPR is the number of falsely detected links divided by the number of all possible false links and 100. The summary of the experiments, i.e. the overall false positive rate (FPR) and the distributions of the link's true positive rate (TPR) across sites are given in Figs. 2 and 3, respectively. The blue violin plots always report the case of normal distributed (non-heteroscedastic) intrinsic noise and the corresponding orange violin plots summarise the case of heteroscedastic noise. The effect of heteroscedasticity and seasonality is then assessed by comparing the distributions obtained from the baseline dataset to the results of the seasonality dataset.

Figure 1 The artificial datasets are generated with a prescribed interaction structure (true network), which is obtained by fitting the test model to the FLUXNET sites. Here we show for four time series lengths the process graphs estimated via both lagged correlation and PCMCI. The data used stems from the homoscedastic realisation of the seasonality dataset of the Hainich site. The significance level was set to 0.01. The number of time lag labels were limited to five in the correlation networks. However, for the longest time series typically the whole range of lags (0–25) was significant.

Figure 2 The distribution of false positive detection rates estimated for the baseline dataset, the seasonality dataset and the anomalised seasonality dataset (mean seasonal cycle subtracted). The distributions are given for different time series length (number of datapoints). Additionally, the distributions are split to show the impact of heteroscedastic noise (orange) compared to normal distributed noise (blue). The significance level of 0.01 is given by a blue horizontal line.

Figure 3 The distribution of the true positive detection rates for each link in the test model estimated for the baseline dataset, the seasonality dataset, and the anomalised seasonality dataset. The distributions are given for different time series lengths (number of datapoints). Additionally, the distributions are split to show the impact of heteroscedastic noise (orange) compared to normally distributed noise (blue).

The FPR of homoscedastic time series in the baseline dataset is in the expected range of 0.01, the chosen significance level, indicating a well-calibrated test due to fulfilled assumptions. The assumption of stationarity is violated as soon as heteroscedasticity or seasonality is present. The effect on the FPR is an increase above 0.01 for time series length of 1 and 5 years, with a much stronger increase due to seasonality (factors of 4 and 8, respectively) than for heteroscedasticity (factor of 2).

The effect of non-stationarities on the TPR differs among the links. The detection of linear links (ℛ_g→𝒯 and 𝒯→𝒯) is not affected by seasonality and slightly improves for heteroscedastic dynamical noise. The detection of non-linear links is improved by seasonality, with the strongest effects in the link 𝒯→𝒢𝒫𝒫. The link 𝒯→ℛ_eco has a stronger non-linearity, and therefore the detection rate shows a weaker effect on seasonality. Furthermore, the coupling coefficient c₅, the base of T, can be close to or be exactly 1 (see Fig. A3). This would actually cause, on the one hand, 𝒯→ℛ_eco to vanish, rendering a detection impossible and, on the other hand, result in a linear dependence of 𝒢𝒫𝒫 on ℛ_eco, which improves its detection. Heteroscedasticity seems to have a slight negative effect on non-linear links. In general, the TPRs in the seasonality dataset are quite high, even for non-linear links, and predominantly above 80 % and often reaching 100 %.

Comparing the TPRs of the non-linear links shows some disparity. The links 𝒯→𝒢𝒫𝒫 and 𝒯→ℛ_eco experience zero detection in the baseline dataset but partially considerable rates in the seasonality dataset with a strong dependence on the time series length. On the contrary, the median of the TPR of the links ℛ_g→𝒢𝒫𝒫 and 𝒢𝒫𝒫→ℛ_eco is above 95 % in the seasonality dataset, even for time series length as short as 91 d, but remains high in the baseline dataset. The removal of the seasonal cycle keeps the TPRs largely unaffected but reduces the FPR. Nevertheless, it still remains above the significance level by factors of 4 and 2 for 5-year and 1-year time series lengths, respectively.

In summary, the seasonality dataset exhibits high TPR even for non-linear links. Compared to stationary time series, the detection of non-linear links actually benefits from seasonality. The high detection, though, comes at the cost of a high false positive rate for time series length of and above 1 year. To a certain degree, the increase in FPR can be counteracted by anomalisation.

3.2 Majadas de Tiètar dataset

First, we look at the link consistency by comparing networks that were obtained for each tower within a month. The comparison is done for 2 months with strongly differing climate conditions: April and August. In Fig. 4 we compare the estimated link strengths (effect size estimated via partial correlation) as long as the corresponding links are significant in at least one network. The confidence intervals are overlapping for the majority of links, suggesting that the uncertainty of the fluxes is much smaller than the observed effects (El-Madany et al., 2018). Exceptions are found for only a few links ($R_{\mathrm{g}}\stackrel{\mathrm{0}}{⟶}T$, $R_{\mathrm{g}}\stackrel{\mathrm{1}}{⟶}LE$, $\mathrm{VPD}\stackrel{\mathrm{1}}{⟶}\mathrm{VPD}$, $H\stackrel{\mathrm{2}}{⟶}H$, $\mathrm{NEE}\stackrel{\mathrm{0}}{⟶}LE$, $H\stackrel{\mathrm{2}}{⟶}\mathrm{NEE}$, $\mathrm{NEE}\stackrel{\mathrm{4}}{⟶}\mathrm{NEE}$; the number above the arrow indicates the lag) where the detection rates do not or barely overlap. Cross-links (a link from one variable to another) with two or more significant appearances are predominantly at zero lag. Approximately half of the links with lag 1 are auto-dependent links (a link from the past of a variable to its present). Comparing the links between the months of April and August, distinct differences can be noticed. First, August has slightly fewer significant links compared to April. Second, the only links remaining that are significant in two or three towers are between atmospheric variables. Third, the remaining link strengths tend to be weaker in August than in April.

Figure 4 Comparison of the networks of three eddy covariance measurement stations (LMa, LM1, LM2) located in Majadas de Tiètar (Spain). Links that are found to be significant in one of the three networks are included. For each link, the calculated strength of all three networks is plotted with its 90 % confidence interval. The colours blue, orange, and green correspond to the towers LMa, LM1, and LM2, respectively. The significance threshold is 0.01. If a link does not pass the significance, it is marked by a black dot. The links are grouped into lag 0 (top), lag 1 (middle), and all lags greater than 1 (bottom). Negative NEE is associated with carbon uptake by the ecosystem. Links at lag 0 are left undirected (–), yet as R_g is set as main driver, links incorporating R_g at lag 0 are directed (→).

The difference among the seasons is further investigated in Fig. 5, which shows process graphs for each month of the year. We combine the networks of the three towers to one process graph by plotting every link that is significant in at least one tower. The process graphs in Fig. 5 visualise clearly gradual changes within the interaction structure of the biosphere–atmosphere system during the course of a year. During the main growing season from February to May, NEE is coupled strongly to the energy fluxes latent (LE) and sensible heat (H). These connections weaken, disappear, or even switch sign during the dry season. Less regularly, NEE also shows connections to radiation (R_g) and temperature (T). Between the atmospheric variables, a basic network between VPD, T, R_g, and H remains intact and relatively constant in strength. The dominance of contemporaneous links is found as well, as already seen in Fig. 4. Aside from the decoupling of NEE from any variable in the dry period, there are additional interesting patterns. For example, the positive reappearance of the link between NEE and LE in September. Here, the onset of precipitation events (see Fig. A10) occurred, which led to strong respiration peaks (Ma et al., 2012). Creating such a network via lagged correlation would result in much more significant links (causing the network to be not interpretable as opposed to PCMCI) and NEE does not decouple from the atmosphere in August (see Fig. A7).

Figure 5 To visualise the gradual changes in interaction structure the networks of the three towers are combined for each month. The number of significant occurrences of a link is given by its width. The link strength, given by the link colour, is calculated by averaging the significant links of the towers. The link's lag is shown in the centre of each arrow, sorted in descending order of link strength. The resulting graphs are shown for April 2014 until March 2015. The significance threshold is 0.01. The networks of April and August, illustrated in Fig. 4, are highlighted by a box.

The above results demonstrate that PCMCI is sensitive enough to capture seasonal differences and certain physiological reasonable biosphere behaviour. Moreover, PCMCI yields a better interpretable network structure than pure correlation approaches.

3.3 Gridded global dataset

The significant lags and MCI values of each climatic variable on NDVI were subject to inspection. In Fig. 6 the maximal MCI value and the corresponding lag are plotted for the links $X_{t-\mathit{\tau }}\to \mathrm{NDVI}:\mathit{\tau }\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$, with X being one of the climatic drivers: radiation, temperature, or precipitation. The chosen significance threshold was set to 0.05. Figure 7 shows the climatic driver with the largest MCI per grid point. PCMCI detects a regionally varying influence of climatic drivers. As expected, the boreal regions are strongly driven by temperature instantaneously, while (semi-)arid regions, which correspond predominantly to grass- or prairie-dominated areas, respond strongest to precipitation at a time lag of 1 month. Radiation is found to have a comparatively low spatial effect, with hotspots in southern and eastern China, central Russia, and eastern Canada.

Figure 6 Influence of climatic drivers on NDVI as calculated by PCMCI. The first column shows the estimated causal influence given as maximal absolute MCI value of climatic drivers on NDVI. The second column gives the time lag at which the maximal absolute MCI value occurs (in months).

Figure 7 Map of the strongest climatic driver (largest absolute MCI value) per grid point.

The dominant lags are found to be 0 and 1. Only a very small fraction of the total area shows a maximal MCI value at a higher lag of two or three months. The lags are also not equally distributed among the climatic drivers. Radiation and temperature are predominantly strongest at lag 0, while precipitation has a much larger fraction of area showing the strongest response at lag 1. Regions where the impact of R_g on NDVI is strongest at lag 1 tend to respond negatively to R_g but positively to precipitation at lag 1. On the other hand, a large part of the regions with the strongest impact of precipitation at lag 0 respond negatively to it but positively to radiation.

In summary, PCMCI estimates coherent interaction patterns that match well with anticipated behaviour based on vegetation type and prevailing climatic conditions.

4 Discussion

Causal discovery methods promise an improved understanding and can help to come up with new hypotheses about the interaction between biosphere and atmosphere (Christiansen and Peters, 2020; Runge et al., 2019a). However, the underlying assumptions need to be properly taken into account. The coupled biosphere–atmosphere system possesses several challenges that potentially violate the underlying assumptions of causal discovery in general and the employed method's assumptions in particular. Here, we investigate the effect of a violation of assumptions on PCMCI network estimates.

With regard to expected non-linearities in biosphere–atmosphere interactions, using a linear independence test within the PCMCI framework may not be adequate. We justify our choice with the following arguments: first, non-linearities are often approximated linearly. Second, a linear regression based test has a much higher power for detecting linear links than a non-parametric test (Runge, 2018) and can hence already detect links at smaller sample sizes. Third, linear partial correlation is easily interpretable, for example, positive and negative MCI values. This motivation is supported on the one hand by the results of the test model and on the other hand by additionally performed analyses on the observational datasets using Gaussian regression and distance correlation as an independence test. These results (see Figs. A8, A9, A11) show similar patterns, but due to the low sample sizes exhibit lower statistical significance. In general, the derived results show a high detection power with a strong consistency in calculated effect strengths on eddy covariance data and global, regularly gridded reanalysis data that leads to easily interpretable patterns. Observed drawbacks are a high FPR in cases of violated assumptions, especially strong periodicity, as well as the appearance of contemporaneous lags in measurement datasets.

4.1 Lessons learned from the test model

The probability of detecting a link with PCMCI depends strongly on a link's MCI effect size, which is larger for strong variance in the driver and a low variance in the receiver (see Sect. 2.5). Several results can be explained by this observation. First, the variance of three out of five drivers of cross dependencies in the test model are either directly or indirectly (via 𝒢𝒫𝒫) influenced by ℛ_g, which has the highest variance of all variables. Consequently, the detection power of the three links is large, almost 100 %. In comparison, the other variables' variances are weaker, since they are influenced by 𝒯, which results in a lower detection power. This is the origin of the disparity in detection rates of the non-linear links. Second, the partially strong increase in TPR of non-linear links (influenced in a multiplicative way by 𝒯) from the baseline dataset to the seasonality dataset can also be explained by this increase in variance. A multiplicative link is actually not generally expected to be found by ParCorr (Runge, 2018), but the value of the multiplicative factor is dominated by the seasonal value, and not the dynamical noise, which might cause a scaling of the dynamical noise terms rather than a random distortion. Third, the dependence on the variance ratio can also explain the difference in TPR between homoscedastic (equal error variance) and heteroscedastic (error variance changing over time) time series, i.e. the variance of ℛ_g and 𝒢𝒫𝒫 exhibits a strong seasonality, with its peak in summer, while the variance of 𝒯 is rather constant. This explains, for example, the strong decrease in TPR for the link 𝒯→𝒢𝒫𝒫 at 91 d time series length when comparing homoscedasticity to heteroscedasticity. The decrease in TPR is less pronounced when another season, implying a different variance, is chosen for this comparison. As links with weak driver variance and strong response variance are more likely to be missed, one may ask what effect this will have on the detection of feedback loops where one variable has low and the other high variance. Here lies a limitation of the test model where no feedback loops were implemented.

Seasonality and heteroscedasticity constitute violations of the stationarity assumption underlying the independence test ParCorr. Seasonality constitutes a common driver in this model. In general, such common drivers increase the dependence among the variables and hence lead to a higher detection rate for true links (TPR) as well as a higher false positive rate (FPR) for absent links if this driver is not conditioned out properly. This additionally causes the TPR and the FPR rate to increase in the seasonality model. As shown in Runge (2018), including the cause of the non-stationarity as an exogenous driver in the analysis allows PCMCI to regress out its influence on the other variables. However, for ParCorr this is only valid if the dependence on the non-stationary driver is linear. Therefore, the regression on ℛ_g fails for 𝒢𝒫𝒫 and ℛ_eco in the test model. With this ill-posed setting, the probability to detect false links increases with increasing time series length or when more periods are included. Stationarity in the mean is obviously also not fully guaranteed when subtracting the seasonal mean. Here we observe that the FPR stays above the significance level for the anomalised seasonality dataset. One can ask whether the FPR stays above the significance threshold because subtracting the seasonal mean does not remove the heteroscedasticity. However, we attribute this high FPR to a not fully removed seasonality since the FPR of both homoscedastic and heteroscedastic time series decreases by roughly the same amount in the anomalised seasonality dataset and the effect of heteroscedasticity is rather weak in the baseline dataset. The increasing FPR with increasing time series length can further raise doubts regarding the analysis of long time series. For such an analysis, though, the assumption of causal stationarity should first be assessed. For example, the link from radiation to GPP vanishes in winter as there is mostly no active plant material left. To account for causal stationarity, the analysis should be limited to time series sections where the causal structure is expected to be similar. This is typically done by limiting the analysis to a specific time period (i.e. “masking”), e.g. a specific season, month, or time of the day. Such masking additionally further reduces influences of remaining seasonality or heteroscedasticity. One can argue, as in Peters et al. (2017), that the causality of a system is invariant even between seasons because the physical mechanism is the same in all seasons. Yet, while the physical, i.e. functional relationship might be constant over time, its imprint in the time series might vary. For example, a functional dependence f(x) might be “flat” for small values of x and linearly increasing for larger values. If only small values occur in the winter season, then the link is absent, while it “appears” only in the summer season. Across all seasons, this can be considered a non-linear functional dependence f(x). In practice, restricting an analysis to different seasons can help in interpreting the mechanism and is performed here using a linear framework.

Summarising the results of the test model, the different detection rates, disparity among non-linear links, and the detection of multiplicative links are largely explainable via the effect of the variance on the link detection. Yet, the discussion revealed the need for further research in several aspects. On the one hand, feedback loops are not included in the test model yet are an important aspect in natural systems. On the other hand, removing non-stationarities is essential to keep the false positive rate in the expected range, but standard procedures of subtracting the mean seasonal cycle are not sufficient. Further, the effect of non-stationarity on the causal network structure needs to be investigated.

4.2 Causal interpretation of estimated networks from observational data

In both the half-hourly time-resolved eddy covariance data and the monthly global dataset, the predominant type of dependence found is contemporaneous. PCMCI leaves these undirected since no time order indicating the flow of causal information is available. Further, as discussed in Sect. 2.2, contemporaneous common drivers or mediators are not accounted for. The consequence is that both spurious contemporaneous and spurious lagged links can appear if they are due to contemporaneous variables. For interactions that are contemporaneous in nature, since they occur on considerably shorter timescales than the time resolution, PCMCI is not the optimal choice regarding a causal interpretation and other methods should be applied in conjunction (Runge et al., 2019a). Further, we faced a trade-off between fulfilling causal assumptions and detection power. In practice, accounting for causal stationarity (by limiting the analysis to certain periods of the dataset) means decreasing the number of available data points but accounting for causal sufficiency leads to an increase in dimensionality by adding variables and increasing the maximal lag. Both will lead to a decrease in detection power, which can affect the network structure. PCMCI alleviates the curse of dimensionality by applying a condition selection step, but one still cannot indefinitely add more variables. Another important factor that affects detection power and dimensionality is the time resolution. There are several points in favour and against increasing time resolution. On the one hand, increasing time resolution can resolve contemporaneous links and potentially increases the detection power due to an increased number of datapoints. On the other hand, the dimensionality increases if the maximal lag is adapted. Further, causal information might be split apart and distributed over more lags, rendering the links at each individual lag less detectable. This can cause links to disappear, but links can also appear if new processes are resolved at a higher timescale. Finally, observational noise (measurement errors) might be larger in higher-resolution data than in aggregated data, as it is averaged out in the latter and thus affects link detection less. Consequently, when comparing network structures based on different settings, i.e. maximal lag, included variables, time resolution, and considered time period, the (dis)appearance of single links among specific variables can stem from several factors, i.e. a change in detection power, a changed (conditional) dependency, or due to a common driver. These factors together with a non-zero false positive detection rate are challenging for a causal interpretation. Therefore, detected links should be interpreted with care and can give rise to new hypotheses and analyses involving further variables. Generally, a causal interpretation is more robust regarding the absence of links (see Sect. 2.2). Specifically, it does not require that all common drivers are observed.

Nevertheless, robust patterns were identified in our studies that are also consistent with other studies. Furthermore, a causal analysis has the advantage of an enhanced interpretability compared to correlative approaches. First of all, we could show that the networks' estimated link strengths are consistent for observational data, even though measurement error affects the data. The dataset used was suitable for this analysis, as the measurement stations are located in a reasonably homogeneous ecosystem that shows only little spatial variation (El-Madany et al., 2018). Thus, also the interaction between biosphere and atmosphere is expected to change only marginally across space within this ecosystem. Second, the gradual changes in plant activity that are taking place in the ecosystem of Majadas de Tiètar throughout the year do emerge very well in the coupling strength of daytime NEE to the atmospheric variables. The observed decoupling during the dry season is in accordance with the one of a soybean field during drought conditions observed by Ruddell and Kumar (2009). The gradual changes in ecosystem activity are not visible in a pure (lagged) correlation analysis or are only visible in colour or density changes but the large number of significant links prevents any detailed interpretation on the physical mechanisms and changes thereof. The large number of significant links compared to the PCMCI networks stems solely from the absence of conditioning in common drivers or mediating variables, which often leads to an overestimation of the link strength in correlation networks. As a result, processes, such as the decoupling of NEE during the dry period, stay hidden. To reduce the effect of confounding, analyses often utilise partial correlation (see, e.g. Buermann et al., 2018). However, a partial correlation can introduce new dependencies (as opposed to removing them) if one conditions on causal effects of the variables under consideration (the “marrying parents” effect). This issue is avoided in PCMCI by only conditioning on past variables. Additionally, PCMCI chooses only relevant variables as conditions by applying the PC condition selection step, which is especially valuable in high-dimensional study cases and improves detection power and computation time (Runge, 2018).

The global study of climatic drivers of vegetation shows a general pattern of lags and dependence strengths of vegetation on climatic variables that is easily interpretable. The boreal regions appear energy limited and especially driven by temperature (see Fig. 6c), while the strongest dependence of (semi-)arid regions on precipitation reflects their limitation in water supply. Two recent studies performed a similar analysis. Both Wu et al. (2015) and Papagiannopoulou et al. (2017b) investigated lagged effects and dependence strengths of NDVI on precipitation, temperature, and radiation. Wu et al. (2015) estimated the lags of the strongest effects via an univariate regression of the climatic drivers on NDVI and subsequently used those lags to fit a multivariate regression model of the climatic drivers on NDVI and determined their relative effects. Papagiannopoulou et al. (2017b) applied a non-linear Granger causality framework utilising a random forest predictive model; the method was presented by Papagiannopoulou et al. (2017a). We recognise that similar patterns are observed in Wu et al. (2015), but the lags at the maximal MCI value are usually lower than the one found in Wu et al. (2015), which stems from the methodical differences. Besides having used anomaly values, PCMCI regresses both NDVI and the climatic drivers on their parents before calculating the MCI value (see Sect. 2.5). This especially removes the influence of autocorrelation. Runge et al. (2014) shows how autocorrelation affects the correlative lag causing it to be larger for stronger autocorrelation; thereby, the correlative lag may become larger than the causal lag. Therefore, according to Fig. 6b, the causal information embedded in monthly resolution is predominantly received within 1 month. Finding the strongest causal links at a time lag up to 1 month appear in agreement with Papagiannopoulou et al. (2017b). Additionally, the spatial distribution of the strongest climatic influences compares well, but there are certain noteworthy differences that do not necessarily stem from masking differences, i.e. that we took only values belonging to the growing season while Papagiannopoulou et al. (2017b) took the whole time series. First, there is little significant Granger causality of water availability found in boreal regions while there are significant negative causal dependencies detected via PCMCI. Second, NDVI in arid regions is not (or is only partially) Granger caused by radiation and temperature but in part shows a negative PCMCI value for those variables. There might be physiological reasons that can explain the PCMCI patterns, i.e. waterlogging or temperatures that are too high. To explain the differences though, we could identify two possible reasons. First, Papagiannopoulou et al. (2017b) masked out negative influences of radiation, arguing that radiation is not negatively affecting NDVI. They found that negative influences of R_g are usually a consequence of poor conditioning on other variables. Second, a precipitation event in boreal regions coincides with a reduction in radiation and temperature. Boreal regions usually do not suffer from water shortages. Thus, they respond more strongly to the reduction of radiation and temperature than precipitation. As precipitation is coupled negatively to radiation and temperature at lag 0, the effect of precipitation on NDVI is found to be negative. Thus, the link $P\stackrel{-}{\to }\mathrm{NDVI}$ might be an effect of the contemporaneous common driver scheme $P\stackrel{-}{\to }{R}_{\mathrm{g}}\stackrel{+}{\to }\mathrm{NDVI}$ and therefore would not be causal. In fact, a similar argumentation can be given for the negative impact of temperature and radiation on NDVI in arid regions.

In summary, we pointed out the need for careful interpretations in applying causal discovery methods and especially highlighted the challenges linked to the study of biosphere–atmosphere interaction via PCMCI. We demonstrated that the network structures estimated from observational data are explainable with respect to plant physiology and climatic effects. Finally, our study shows that causal methods can deliver better interpretability and a much improved process understanding in comparison to correlation and bivariate Granger causality analyses that are ambiguous to interpret since they do not account for common drivers.

4.3 Outlook

The preceding discussion has shed light on the merits of PCMCI, as well as the challenges of applying causal discovery methods. Runge et al. (2019a) discuss further challenges and methods and give an outlook of how multiple methods can be combined to alleviate limitations.

5 Conclusions

Here we tested PCMCI, an algorithm that estimates causal graphs from empirical time series. We specifically explored two types of datasets that are highly relevant in biogeosciences: eddy covariance measurements of land–atmosphere fluxes and global satellite remote sensing of vegetation greenness. The causal graphs estimated from the eddy covariance data collected in a Mediterranean site confirm patterns we would expect in these ecosystems: during the dry season's plants senescence, for instance, the ecosystem's carbon cycle (NEE) decouples from meteorological variability. On the contrary, during the main growing season, with warm and humid conditions, strong links between NEE, LE and H characterise the graph. Using the causal framework, not only the strongly contrasting states emerge in the graph structure but also the gradual transitions that relate to minor changes like the connectivity of sensible heat to temperature with progressing dryness. A purely correlative analysis, instead, is not able to resolve these patterns. PCMCI allows us to identify and focus on fewer (only those that are highly relevant) dependencies. Applying the approach to three replicated eddy covariance systems shows the robustness of the method to random errors in the fluxes measurements and confirms one of the assumptions of eddy covariance: above a relatively homogeneous terrain the fluxes measured should be spatially invariant as should the underlying causal relationship between climate and fluxes. The global analysis of NDVI in relation to climatic drivers confirms the known patterns of dependence strengths of vegetation on climatic variables: boreal regions are energy limited and especially driven by temperature and secondarily by radiation, while in semi-arid regions vegetation dynamics are strongly dependent on water supply. However, obtained response times of vegetation to climatic variations are lower using PCMCI than correlation which can be attributed to a better treatment of the autocorrelation in the time series and cross-relations among climate variables. Compared to merely correlative approaches, this leads to a interpretable pattern of driver–response relationships. In short, the new developments achieved in causal inference allow us to gain well-constrained insights on processes that would otherwise be drowning in the correlation chaos. Therefore, we hope that this study fosters usage of causal inference in analysing interactions and feedbacks of the biosphere–atmosphere system and furthermore exhibits our demand of further developments.

Appendix A

Table A1 PCMCI parameters that were used differently from default settings.

Wohlfahrt et al. (2008)Lund et al. (2012)Meyer et al. (2015)Suni et al. (2003)Beringer et al. (2011a)Thum et al. (2007)Hutley et al. (2011)Delpierre et al. (2016)Cernusak et al. (2011)Berbigier et al. (2001)Beringer et al. (2007)Rambal et al. (2004)Beringer et al. (2011b)Bonal et al. (2008)Leuning et al. (2005)Vitale et al. (2016)Moureaux et al. (2006)Valentini et al. (1996)Aubinet et al. (2001)Marcolla et al. (2003)Saleska et al. (2003)Marcolla et al. (2011)Brooks et al. (1997)Spano et al. (2004–2014)Bond-Lamberty et al. (2004)Rey et al. (2002)Wang et al. (2002a)Chiesi et al. (2005)Wang et al. (2002b)Galvagno et al. (2013)Wang et al. (2002c)Matsumoto et al. (2008)Chen et al. (2006)Jacobs et al. (2007)Rayment and Jarvis (1999a)Kurbatova et al. (2008)Rayment and Jarvis (1999b)Fischer et al. (2007)Merbold et al. (2014)Schade et al. (1999)Zielis et al. (2014)Wofsy et al. (1993)Imer et al. (2013)McDowell et al. (2004)Etzold et al. (2011)Ruehr et al. (2012a)Ammann et al. (2009)Pryor et al. (1999)Dietiker et al. (2010)Gitelson et al. (2003)Dušek et al. (2012)Cassman et al. (2003a)Bernhofer et al. (2009–2014)Cassman et al. (2003b)Anthoni et al. (2004)Ruehr et al. (2012b)Prescher et al. (2010a)Scott et al. (2008)Knohl et al. (2003a)Tang et al. (2003)Prescher et al. (2010b)Hatala et al. (2012)Lindauer et al. (2014)Rothstein et al. (2000)Bernhofer et al. (2008–2014)Nave et al. (2011)Bernhofer et al. (2010–2014)Xu et al. (2004)Grünwald and Bernhofer (2007)Scott et al. (2006)Westergaard-Nielsen et al. (2013)Emmerich (2003)Pilegaard et al. (2011)

Table A2 List of FLUXNET sites used for the generation of artificial datasets and the time period used.

Scheme A1 Schematic overview of the aim and approach of PCMCI using time series graphs as a visualisation tool.

Figure A1 Distribution of coupling coefficients obtained after fitting the test model to the Fluxnet sites. Shown here are the distributions used for generation of heteroscedastic time series.

Figure A2 Distribution of coupling coefficients obtained after fitting the test model to the Fluxnet sites. Shown here are the distributions used for generation of homoscedastic time series.

Figure A3 Distribution of time lags obtained after fitting the test model to the Fluxnet sites. Shown here are the distributions used for generation of heteroscedastic time series.

Figure A4 Distribution of time lags obtained after fitting the test model to the Fluxnet sites. Shown here are the distributions used for generation of homoscedastic time series.

Figure A5 Observed (blue) and test model (orange) time series for the Hainich Fluxnet site. The model data were produced with heteroscedastic noise.

Figure A6 The same as Fig. 5 but using simple correlation analysis to estimate the graph structures. The number of significant occurrences of a link is given by its width. The link strength, given by the link colour, is calculated by averaging the significant links of the towers. Link labels indicating the lag were removed to improve link visibility. They typically ranged from 1 to 8 (full range of possible lags). The resulting graphs are shown for April 2014 until March 2015. The significance threshold is 0.01.

Figure A7 The same as Fig. 4, but the analysis was performed using a non-linear independence test. Comparison of the networks of three eddy covariance measurement stations (LMa, LM1, LM2) located in Majadas de Tiètar (Spain). Links that are found to be significant in one of the three networks are included. For each link, the calculated strength of all three networks is plotted with its 90 % confidence interval. The colours blue, orange, and green correspond to the towers LMa, LM1, and LM2, respectively. The significance threshold is 0.01. If a link does not pass the significance, it is marked by a black dot. The links are grouped into lag 0 (a, b), lag 1 (c, d) and all lags greater than 1 (e, f). Links at lag 0 are left undirected (−), yet as R_g is set as main driver, links incorporating R_g at lag 0 are directed (→). Note that Gaussian process regression and distance correlation (GPDC) only yields positive link strengths. Further, the strength values estimated with GPDC are rather weak due to the low number of datapoints and the larger sensitivity of that method to the sample size.

Figure A8 The same as Fig. 5, but the analysis was performed using a non-linear independence test. The number of significant occurrences of a link is given by its width. The link strength, given by the link colour, is calculated by averaging the significant links of the towers. The link's lag is shown in the centre of each arrow, sorted in descending order of link strength. The resulting graphs are shown for April 2014 till March 2015. The significance threshold is 0.01. Note that GPDC only yields positive link strengths. Further, the strength values estimated with GPDC are rather weak due to the low number of datapoints and the larger sensitivity of that method to the sample size

Figure A9 Daily aggregated precipitation in Majadas de Tiètar measured at the three tower sites from April 2014 to March 2015. Missing values are plotted as gaps.

Figure A10 Influence of climatic drivers on NDVI as calculated by PCMCI in conjunction with the non-linear independence test GPDC (similar to Fig. 6). The first and second columns show the estimated causal influences of climatic drivers on NDVI at lag 0 and 1, respectively.