Inferring causal associations in hydrological systems: A comparison of methods

Inferring causal associations in hydrological systems: A comparison of methods | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Inferring causal associations in hydrological systems: A comparison of methods Hanxu Liang, Wensheng Wang, Bin Chen, Li Guo, Hu Liu, Siyi Yu, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4643196/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 16 Apr, 2025 Read the published version in Stochastic Environmental Research and Risk Assessment → Version 1 posted 9 You are reading this latest preprint version Abstract Many research issues in hydrological systems are intrinsically causal, aiming to determine whether and how one factor affects another. Although causal inference methods have been applied more or less in hydrology, there still remains a lack of systematic comparison between different methods. Here, four popular methods in the causal inference community, including the cross-correlation function (CCF), convergent cross mapping (CCM), transfer entropy (TE), and a causal network learning algorithm (PCMCI+) were selected, with a detailed explanation of their basic principles and underlying assumptions. Next, the performances of these methods were evaluated in large sample tests and sensitivity analysis using synthetic time series generated by a conceptual hydrological model with two predesigned causal structures. Then, the four methods were applied in two real-world cases to further understand their characteristics. The findings show the superior performance of the PCMCI + method in synthetic cases and a commendable level of interpretability in real cases, thus warranting its broader application in hydrological systems. The limitations of the other three methods, especially in effectively addressing confounding and mediating factors, led to several unreasonable causal links. Furthermore, the emergence of conflicting results among different methods in real-world applications underscores the necessity for a multifaceted understanding based on their particular assumptions and constraints. A comprehensive application of diverse methods according to the specific issue is encouraged for the robustness of conclusions, with their assumptions clearly stated in advance. Overall, our research reveals the potential and limitations of different causal inference methods in comprehension of complex interactions within hydrological systems, serving as a useful guide for their further prosperity in hydrology. causal inference cross-correlation function (CCF) convergent cross mapping (CCM) transfer entropy (TE) PCMCI+ hydrological systems Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Figure 12 1. Introduction The hydrological systems, which encompass the interaction of multiple variables (precipitation, runoff, evaporation, etc.) across spatiotemporal scales, are complex giant systems with both deterministic and stochastic dynamics (Sang et al., 2015 ). Many hydrological questions are inherently causal, which aim to understand whether and how the cause variables impact the effect variables. For instance, how to identify and quantify interactions among variables in the hydrologic cycle (Good et al., 2015 ; Kleidon and Renner, 2013 ); how to select appropriate influence factors for robust hydrologic predictions (Apaydin and Sibtain, 2021 ; Wang et al., 2023 ; Chen et al., 2011 ); how to determine the direct and indirect factors in the analysis of streamflow or hydrologic extremes (i.e., flood and drought) under the changing environment (Liang et al., 2023 ; Zhang et al., 2021a ; Peng et al., 2024 ). It is widely acknowledged that correlation does not imply causation (Altman and Krzywinski, 2015 ). Associations can arise between variables in the presence (i.e., X causes Y ) and absence (i.e., they have a common cause) of a causal relationship, as stated in Reichenbach’s common cause principle (Reichenbach, 1956 ). Relying solely on correlation-based methods may not fully encapsulate the intrinsic causal mechanisms in hydrological systems, especially within the complex spatiotemporal fabric characterized by numerous variables. Besides, Blöschl et al. ( 2019 ) recently stated in the 23 unsolved problems for hydrologic sciences that“ Questions remain focused on process-based understanding of hydrological variability and causality at all space and time scales ”, which highlights the urgent need for the causal perspective and framework to promote a better understanding of complex hydrological systems. Generally, replicated interventional experiments on variables provide an explicit way to understand causal processes. However, it is often unfeasible or unethical in the real world, such as operating and intervening in regional or global hydrological variables. Additionally, with the development of monitoring technology, several experimental catchments have been established in recent years (Zhang et al., 2021b ). However, the construction and maintenance are costly, and due to the spatial heterogeneity and scale problems of hydrological systems (Bergström and Graham, 1998 ; Blöschl et al., 2019 ), the extensibility and upscaling of discovered mechanisms are also challenging. Moreover, data-informed computer simulations may provide alternative randomized controlled experiments, while these are often time-consuming and inaccurate, and demand massive expert knowledge, which in turn may impose strong mechanistic assumptions on the system (Runge et al., 2019a ). Fortunately, with the development of observational techniques and data sciences, the amount of available time series data in hydrologic sciences has been increasing over the years, which comes from satellite remote sensing observations, station-based or field site measurements, earth system modeling, and reanalysis products (Li et al., 2023a ). Such data repositories, along with growing computational efficiency, open up an alternative avenue to use data-driven methods without the intervention of hydrological systems, namely observational causal inference. Over the past decades, with the progress of statistics and computer science, methods for causal inference from observational time series have been developed rapidly. The most classic work can be traced back to the Granger causality (GC; Granger, 1969 ), which approves the causal impact from the variable X to Y if the past of X contributes to the prediction of the future of Y. However, the validity of the GC framework remains controversial since it only reveals the statistical associations via linear vector autoregression (Shojaie and Fox, 2022 ). Later, the nonlinear extension of the original GC framework was developed by integrating information theory, that is, the transfer entropy (TE; Schreiber, 2000 ). Except for the GC framework, some complementary perspectives emerged, from the nonlinear dynamics system based on the reconstruction of attractors, known as the convergent cross mapping (CCM; (Sugihara et al., 2012 )). More recently, the graph-based causal network learning algorithms, such as the Peter–Clark momentary conditional independence (PCMCI) proposed by Runge et al. ( 2019b ), become increasingly popular in the causal community. Unlike the data-driven machine learning methods that primarily concentrate on classification or prediction (Zounemat-Kermani et al., 2021 ), the purpose of all causal methods is to discover and quantify the causal interactions of the underlying system based on the measured time series (Runge et al., 2023 ). Since explainable machine learning has attracted great interest of researchers recently (Angelov et al., 2021 ), it is still challenging for machine learning, even deep learning models to extract complex mechanisms from this black box (Rudin, 2019 ). Thus, causal inference methods are significant for supplementing predictive machine learning to enhance our theoretical cognition of underlying systems (Reichstein et al., 2019 ). Over the years, causal inference methods have shown huge potential in numerous fields across Earth system science, for example, the regularity discovery, process understanding, hypothesis validation, and physical model improvement (Su et al., 2023 ; (Faybishenko, 2017 ). Unfortunately, in hydrology, traditional correlation and regression methods are still the most commonly used tools (Massei et al., 2006 ; Xu et al., 2022 ; Yu et al., 2011 ), and the application of modern causal methods remains incipient, such as the feedback between rainfall and soil moisture (Luo et al., 2023 ; Wang et al., 2018 ), the interactions between groundwater and streamflow (Bonotto et al., 2022 ; Zhao et al., 2023b ), and the hydrological connectivity in the karstic areas (Delforge et al., 2022 ; Li et al., 2023b ). The modern causal inference methods do have huge potential in hydrological systems if their underlying assumptions and methodological challenges have been fully considered (Runge et al., 2019a ). Since limited exploratory studies on the applications of causal inference methods in hydrology have been conducted recently, there still remains much confusion about how and when to apply them, and how to interpret their results, thus requiring a comprehensive comparison and evaluation between them to guide the further popularization in hydrological sciences. In this study, four popular methods in the causal inference community, namely cross-correlation function (CCF), convergent cross mapping (CCM), transfer entropy (TE), and a causal network learning algorithm (PCMCI+) were selected, with a detailed explanation of their basic principles and underlying assumptions. Then, the performances of the methods were systematically evaluated in both synthetic and real hydrological systems, thus revealing their potential and limitations for application in hydrology. The remainder of this paper is organized as follows. Section 2 provides a comprehensive description of the four causal inference methods in our comparison study. Section 3 successively illustrates the setup of the hydrological model structure, evaluation metrics, and experimental results in the synthetic case, which includes the large sample test and sensitivity analysis. Section 4 presents the performances of the causal inference methods in two real case studies. Further discussion and conclusions are illustrated in Sections 5 and 6 , respectively. 2. Causal inference methods Four causal inference methods, including the CCF, CCM, TE, and PCMCI+, which are popular and representative in the causal inference community (Runge et al., 2019a ; Runge et al., 2023 ), were chosen for our research. The basic characteristics of these methods are listed in Table 1 . Table 1 Summary of the characteristics of the four causal inference methods adopted in this research. Cross-correlation function (CCF) Convergent cross mapping (CCM) Transfer entropy (TE) Causal network learning algorithm (PCMCI+) Linear/nonlinear linear nonlinear nonlinear linear/nonlinear a Type of system stochastic deterministic stochastic stochastic Detection of contemporaneous causal links No Yes/No b No Yes Other assumptions / nonlinear dynamical system / causal sufficiency, causal Markov, causal faithfulness References (Cryer and Chan, 2008 ) (Sugihara et al., 2012 ); (Ye et al., 2015 ) (Schreiber, 2000 ) (Runge, 2020 ) Note: a Depends on the algorithms for conditional independence detection. b Depends on whether the principle of priority of the cause is adopted. 2.1 Cross-correlation function For a driving variable 𝑋 and a response variable 𝑌 , causality can be identified through the cross-correlation function (CCF) with the principle of priority of the cause. As shown in Fig. 1 a, the variable 𝑋 𝑡−d is considered to be the cause of variable Y 𝑡 if the Pearson’s correlation coefficient ρ between X t − d and Y t is significant during their common time domain for at least one value of d up to d max : where Cov ( X t−d , Y t ) represents the cross-covariance at lag time d ; σ x and σ y denote the standard deviations of time series x and y , respectively. The significance of the hypothesis (i.e., 𝜌 is different from zero) is usually evaluated using the Student’s t-test (Student, 1908 ) with a p -value. Given a significance level 𝛼, significant relationships are considered if the p -value is lower than 𝛼. In this study, significant correlations identified by the CCF method are considered as a way to reveal potentially causal links between two time series at a causal delay. However, according to Reichenbach’s common cause principle (Reichenbach, 1956 ), this linear correlation may not imply causation, owing to the interference of confounders, autocorrelation of variables, and nonlinear dynamical associations (Dean and Dunsmuir, 2016 ). Despite these inherent limitations as a bivariate linear method, the CCF method has still been widely applied in many disciplines including hydrology due to its simplicity, efficiency, and linear interpretability (Massei et al., 2006 ; Xu et al., 2022 ; Yu et al., 2011 ). 2.2 Convergent cross mapping Convergent cross mapping (CCM), proposed by Sugihara et al. ( 2012 ), aims to infer causality between two time series in nonlinear dynamical systems using the theory of time delay embedding (Takens, 1981 ). It detects topological similarities between reconstructed attractor manifolds (i.e., consistent behaviors when revisiting similar states). The basic assumption of the CCM is that variables X and Y belong to the same nonlinear dynamical system, and if X causes Y , the cause variable X leaves informative features in the affected variable Y , and thus X can be estimated from its features in Y (Fig. 1 b). Specifically, to identify the causation between variables X and Y , the first stage is to reconstruct the shadow manifold by the lagged time series of each variable. For instance, the trajectories M x ={ X t , X t−τ , X t −2 τ } and M y ={ Y t , Y t−τ , Y t− 2 τ } are the shadow manifold of variables X and Y , respectively, and τ denotes the time delay. Secondly, for a single point in M y , k (= E + 1) nearest neighbors \({Y_{{t_m}}}\) ( m = 1, 2,… k ) are identified based on their Euclidean distance to Y , where E is the embedding dimension. Next, time indices of these identified nearest neighbors are utilized to map corresponding points \({X_{{t_m}}}\) in M x , and their weighted average values considering the Euclidean distances are calculated to estimate X , shown as follows: $${\mathop X\limits^{ \wedge } _t}|{M_y}=\sum\limits_{{m=1}}^{k} {{\omega _m}} {X_{{t_m}}}$$ 2 $${\omega _m}=\frac{{{u_m}}}{{\sum\limits_{{n=1}}^{k} {{u_n}} }}$$ 3 $${u_m}={e^{ - \frac{{d({Y_t},{Y_{{t_m}}})}}{{d({Y_t},{Y_{{t_1}}})}}}}$$ 4 where \({\mathop X\limits^{ \wedge } _t}\) represents the estimated values of X t from M y , and d ( X , Y ) denotes the Euclidean distance between vectors X and Y . Finally, the cross-mapping skill, represented by Pearson's correlation coefficient 𝜌 between observed X t and the estimated \({\mathop X\limits^{ \wedge } _t}\) , quantifies the strength of causality from variable X to Y . Besides, 𝜌 will increase from 0 to a convergent value as the cross-mapping process undergo all points of the shadow manifold seriatim. Referred to Ombadi et al. ( 2020 ), the significance threshold of 𝜌 was set as 0.3. Given that the CCM method is based on the chaos theory, it should be applied under the strict assumptions of the deterministic mathematical models, namely, low-dimensional systems with infinite length, no noise, and acyclic and non-intermittent relations (Khatibi et al., 2012 ; Sivakumar, 2004 ). Additionally, according to Sugihara et al. ( 2017 ), convergence can be achieved when variables are subject to synchrony, a well-known phenomenon resulting from strong coupling or dynamic resonance, which might lead to spurious bidirectional causal relationships, especially for hydrological variables (Ombadi et al., 2020 ; Zhao et al., 2023a ). Therefore, in this study, the time-lagged CCM (Ye et al., 2015 ), which aims to overcome the synchrony with the asymmetric patterns of time dependencies (i.e., the principle of priority of the cause), was adopted for our research. Distinguished from the traditional CCM method that uses Y t to cross map X t , the time-lagged CCM method adopts Y t+d instead. Note that d is the lag time of cross-mapping, ranging from 0 to d max , which is different from the embedding delay τ in reconstruction of the shadow manifold. 2.3 Transfer entropy Transfer Entropy (TE) is an information-theoretic approach that quantifies the temporally asymmetric transfer of information between two time series (Schreiber, 2000 ). It can be deemed as the nonparametric extension of the Granger Causality (Granger, 1969 ) and be adopted to infer the causal association from X to Y if knowing the information about the past of variable X can reduce the uncertainty of variable Y . Due to its no assumptions on the underlying structure of data and applicability to nonlinear separable systems, the TE method has been widely applied in hydrology (Bennett et al., 2019 ; Konapala et al., 2020 ; Moges et al., 2022 ). To calculate TE, the concept of mutual information is introduced (Kinney and Atwal, 2014 ), which can be considered as the overlapping information between two random variables X and Y , shown as follows: $$I(X;Y)=H(X)+H(Y) - H(X,Y)$$ 5 where I(X; Y) denotes the mutual information between variables X and Y ; H is the Shannon entropy (Shannon, 1948 ), a measure of a variable’s average uncertainty; H(X) and H(X, Y) represent the univariate and bivariate information entropy, respectively, calculated as: $$H(X)= - \int\limits_{X} {p(x){{\log }_2}(p(x))dx}$$ 6 $$H(X,Y)= - \int\limits_{Y} {\int\limits_{X} {p(x,y){{\log }_2}(p(x,y))dxdy} }$$ 7 where p represents the probability distribution of the random variable. Then, TE can be understood as the shared mutual information between the past of the cause variable X and the present of the affected variable Y , conditioned on the past of Y (Fig. 1 c), calculated as: $$T{E_{X \to Y}}=I({Y_t};{X_{t - d}}|{Y_{t - 1}})$$ 8 $$I(X;Y|Z)=H(XZ)+H(Y|Z) - H(X,YZ)$$ 9 $$H(X|Y)=H(X) - I(X;Y)$$ 10 where d is the lag time in the transfer of information between X and Y ; I(X;Y|Z) and H(X|Y) denote the conditional mutual information and conditional entropy, respectively. Besides, Y t− 1 was selected to represent the past of variable Y based on the assumption of the Markov process that the immediate history always provides the most information (Budakoti et al., 2021 ; Ruddell and Kumar, 2009 ). To estimate the probability distribution of the variables, following Ruddell and Kumar ( 2009 ), a histogram-based approach was adopted with 11 bins. Considering that the hydrological time series includes several zero values, which may cause deviation in the calculation of entropy, we used the approach proposed by Gong et al. ( 2014 ) to handle the zero and nonzero values separately in the binning. Additionally, the significance of TE was tested by the Monte Carlo methods and Student’s t-test. The TE values were in comparison with a series of 500 randomly shuffled surrogates in which any temporal correlations between the two time series were broken (Moges et al., 2022 ). The TE value is considered significant if it exceeds the threshold, defined as the value corresponding to the 95th percentile of the random samples. 2.4 Causal network learning algorithms The causal network learning algorithms, based on a series of graphical rules that dominate the identification of system causal associations, have been developed for reconstructing large-scale causal graphs with high dimensionality of variables (Runge et al., 2019b ). To infer the structure of the causal graph, three underlying assumptions, namely causal sufficiency, causal Markov, and causal faithfulness are usually required (Runge, 2018 ). The causal sufficiency assumes that there are no other unobserved (or latent) variables that directly or indirectly impact any pair of our variable sets (i.e. all variables are directly observed). The other two assumptions allow to relate d-separations to conditional independence in the graph under this distribution: X d-sep Y | Z ⟺ X ⊥ Y | Z , where causal Markov and causal faithfulness assumptions correspond to the forward and backward arrows, respectively. Two nodes can be deemed as d-separated given Z if and only if all paths are blocked given Z (Pearl, 2009 ). The causal faithfulness assumption precludes the case where X affects Y in two directions with positive and negative effects canceling out. In this study, the latest causal network learning algorithm, namely Peter-Clark momentary conditional independence plus (PCMCI+), which solves the problems for detecting contemporaneous and time-lagged causal links in autocorrelated high dimensional time series (Runge, 2020 ), was adopted. The algorithm can be divided into two stages, namely the skeleton discovery phase and the orientation phase (Fig. 1 d). In the skeleton discovery phase stage, the Peter–Clark (PC) algorithm begins with a fully connected graph and then tests for the elimination of links between variables iteratively by conditioning sets of increasing cardinality (Spirtes and Glymour, 1991). Specifically, it starts with a graph where all unconditionally ( p = 0) dependent variable pairs are connected with the assumption of the stationarity for causal links, and iteratively tests conditional independence with an increasing number of conditions p . Lagged links are oriented forward in time based on the principle of the priority of cause, whereas contemporaneous links are left undirected (circle marks at the ends) in this skeleton discovery phase. For instance, X t −1 and Z t (black nodes) have already been correctly identified as independent in the second iteration step ( p = 1) where the dependence through Y t −1 (blue box) is conditioned out, while we need to condition on two variables to test whether Z t −2 and W t are independent ( p = 2). To further eliminate spurious links for all ordered pairs and increase the detection power, the contemporaneous conditions are iterated over by momentary conditional independence (MCI) tests with partial correlation: \(MCI:X_{{t - d}}^{i} \bot X_{t}^{j}|\mathop P\limits^{ \wedge } (X_{t}^{j})\backslash \{ X_{{t - d}}^{i}\} ,\mathop P\limits^{ \wedge } (X_{{t - d}}^{i})\) (11) where \(\widehat{P}\left({X}_{t}^{j}\right)\) and \(\widehat{P}\left({X}_{t-d}^{i}\right)\) represent the lagged parents of \({X}_{t}^{j}\) and \({X}_{t-d}^{i}\) identified in the previous PC step, respectively; ‘\’ means precluding \({X}_{t-d}^{i}\) from \(\widehat{P}\left({X}_{t}^{j}\right)\) . After the iteration through subsets S \(\subset\) X t of contemporaneous links, the spurious adjacencies can be fully removed. The MCI partial correlation values were used to represent the link strength of causality, with the significance level α pc set to 0.05 in the Student’s t-test. In the second orientation phase, the left contemporaneous links can finally be directed by a series of rules. For example, when finding that W t −1 and Z t are independent conditional on Z t −1 , while not independent conditional on W t , the causal links from Z t to W t can be identified since the other causal directions are not in accordance with the observed conditional independencies. However, there may not exist such rules to distinguish the direction between the X t and Y t due to the Markov equivalence class (Pearl, 2009 ). More details about the PCMCI + algorithm, along with its pseudo code, can be found in Runge ( 2020 ). To the best of our knowledge, the PCMCI + algorithm has not been applied in hydrological systems hitherto, and this study provides the first evaluation of its potential for further applications. 3. Application to synthetic case study 3.1 Construction of conceptual hydrological model To evaluate the ability of the above four methods for identifying causal associations in hydrological systems, firstly, a conceptual hydrological model was employed to generate the synthetic time series, where the underlying causal associations have already been acquired and can be deemed as true relationships. The EXP-HYDRO model (Patil and Stieglitz, 2014 ), which characterizes the complexity of hydrological processes while maintaining as few variables and parameters as possible to ensure the interpretability of the constructed causal network, was chosen as our basic model framework. We simplified and split the original EXP-HYDRO model into two structures, namely the rainfall-runoff structure (Model S1, Fig,2a) and the snowmelt-runoff structure (Model S2, Fig. 2 b). Model S1 includes four variables: effective precipitation P , soil water storage S , interflow I , and runoff Q . Besides, the model contains three parameters: the maximum soil water storage S max , and nonlinear storage-discharge parameters K s , and γ . The model structure and corresponding causal network are shown at the top of Fig. 2 a. The only forcing variable P was generated by a stochastic weather generator (WeaGETS; Chen et al., 2010 ), including a Markov chain for occurrence and a gamma distribution for quantity. Specifically, the occurrence of precipitation was generated with a binary variable \(\widehat {P}\) (i.e., wet or dry days) using the first-order Markov chain (Richardson, 1981 ): \({\widetilde {p}_{i,j}}(t)=probability({\widehat {P}_t}=j|{\widehat {P}_{t - 1}}=i){\text{ ; }}i,j=0,{\text{ }}1{\text{ ; }}t>1\) (12) where \({\widetilde {p}_{i,j}}\) is the transition probability from the state i to j ; the states 1 and 0 represent the wet and dry days, respectively. Considering that precipitation either occurs or not on a given day, \({\widetilde {p}_{0,1}}+{\widetilde {p}_{0,0}}={\widetilde {p}_{1,0}}+{\widetilde {p}_{1,1}}=1\) . In this study, we set \({\widetilde {p}_{0,0}}={\widetilde {p}_{1,1}}=0.8\) , \({\widetilde {p}_{0,1}}={\widetilde {p}_{1,0}}=0.2\) . The values of \({\widetilde {p}_{0,0}}\) and \({\widetilde {p}_{1,1}}\) were set relatively high in order to raise the persistence in the model, which could make some obstacles to the detection of causal associations. Next, the quantity of precipitation was generated by a Gamma distribution: \(\widehat {P}\sim Gamma(\alpha ,\beta)\); Generally, for precipitation data, the parameter α is smaller than 1 while β has a wide range of possible values (Geng et al., 1986 ). In this study, α and β were set as 0.6, and 6, respectively. Besides, adding some noise to the forcing variable helps to approach the real situation and test the anti-interference ability for causal inference methods. Here, \({P_t}={\widehat {P}_t}+{\eta _p}\) , and \({\eta _p}\) is the white Gaussian noise (zero mean and \(\sigma _{P}^{2}\) variance). The \(\sigma _{P}^{2}\) can be calculated as the expectation value of P 2 divided by the signal-to-noise ratio ( SNR ), as follows: \({\sigma _P}^{2}=E{\text{ }}[{P^2}]/SNR\) and the SNR can be transformed into the unit of dB: \(SNR(dB)=10{\log _{10}}SNR\) . The other variables S , I , and Q were determined by the equations of water balance and nonlinear storage-discharge relationship, as: $$\frac{{d{S_t}}}{{dt}}={P_{t - 1}} - {I_{t - 1}}$$ 13 $${I_t}={K_s} \times S_{{t - 1}}^{\gamma }$$ 14 $${Q_t}=\left\{ \begin{gathered} {S_{t - 1}}+{P_{t - 1}} - {S_{\hbox{max} }}{\text{ }};{\text{ }}{S_t} \geqslant {S_{\hbox{max} }} \hfill \\ {\text{ }}0{\text{ ; }}{S_t}<{S_{\hbox{max} }} \hfill \\ \end{gathered} \right.$$ 15 where K s represents the speed of reduction in water storage, and γ is an exponent for S. Model S2 includes four variables: temperature T , snowmelt M , soil water storage S , and interflow I , and three parameters: the nonlinear storage-discharge parameters K s , and γ , and snow-melting parameter K melt . The model structure and corresponding causal network are shown at the top of Fig. 2 b. The only forcing variable T was generated by a Normal distribution with a mean of zero. We also added some noise to this forcing variable: \({T_t}={\widehat {T}_t}+{\eta _T}\) , where \({\eta _T}\) is the white Gaussian noise (zero mean and \(\sigma _{T}^{2}\) variance). The calculation of \(\sigma _{T}^{2}\) is similar to \(\sigma _{P}^{2}\) . The other variables S , I , and M were determined by the equations of water balance and nonlinear storage-discharge relationship, as: $${I_t}={K_s} \times S_{{t - 1}}^{\gamma }$$ 16 $$\frac{{d{S_t}}}{{dt}}={M_{t - 1}} - {I_{t - 1}}$$ 17 $${M_t}=\left\{ \begin{gathered} {K_{melt}} \times {T_{t - 1}}{\text{ }};{\text{ }}{T_{t - 1}} \geqslant 0 \hfill \\ {\text{ }}0{\text{ ; }}{T_{t - 1}}<0 \hfill \\ \end{gathered} \right.$$ 18 The parameters of the conceptual hydrological model are listed in Table 2 and the corresponding frequency distributions of the generated hydrological variables are shown at the bottom of Fig. 2 . The frequency distributions are estimated with a sample size of 2000 and SNR (dB) of 40. Since most hydrological series follow a skewed distribution and have considerable zeros values (Gong et al., 2014 ), especially for the forcing variable rainfall/snow, the synthetic series match these properties and can represent the real hydrological system to some extent. Besides, for each model structure, 100 datasets were generated to ensure the robustness of the assessment. In the sensitivity evaluation of causal methods, sample size varies from 100 to 2000; SNR (dB) from 2 to 40; data missing rate from 10–60% with two scenarios, namely synchronous and asynchronous sparsity. In the former scenario, all variables are missing at the same time indices that are generated randomly, while in the latter, the time indices vary from the variables. The missing values are complemented by the arithmetic mean of observed values of the same variable (Gao et al., 2018 ). Table 2 The parameters of the conceptual hydrological model. Parameter Value maximum soil water storage ( S max ) 20 Model S1 storage-discharge parameter 1 ( K s ) 0.01 storage-discharge parameter 2 ( γ ) 1.5 snow-melting parameter ( K melt ) 2 Model S2 storage-discharge parameter 1 ( K s ) 0.5 storage-discharge parameter 2 ( γ ) 0.6 In this study, three metrics, namely True Positives Rate (TPR), False Positives Rate (FPR), and Accuracy Rate (AR), were adopted to evaluate the performance of causal inference methods: $$TPR=\frac{{TP}}{{TP+FN}}$$ 19 $$FPR=\frac{{FP}}{{FP+TN}}$$ 20 $$AR=\frac{{TP+TN}}{{TP+FP+TN+FN}}$$ 21 where TP and FN denote the true positive and false negative, respectively, which means that the causal link generated by the model is correctly and incorrectly identified, respectively; TN and FP are true negative and false positive, respectively, which means that the non-causal link is correctly and incorrectly identified, respectively. The lower values of FPR and the higher values of TPR and AR indicate the better performance of the tested causal method. 3.2 Performance in large sample tests In this section, 100 datasets were generated with a sample length of 2000 to evaluate the performance of different causal methods. Figure 3 presents the causal structures of the hydrological model identified by four causal methods. The results represent the average behavior of the experiment, i.e., a causal link exists only if it is identified by more than half of the 100 simulations. The detection of the true causal structures in Models S1 and S2 is mainly disturbed by confounding and mediating variables, respectively, which are typical challenges for causal inference. In Model S1, all methods correctly identified the causal links P → Q , P → S , and S → I . The CCF method incorrectly identified the links Q → I and Q → S due to the influence of the confounding variables P (simultaneously affecting S and Q ) and S (simultaneously affecting Q and I ) respectively. Besides, the indirect link P → I results from the mediating variable S . These incorrect links in other causal methods (CCM, TE, and PCMCI+) also share the similar reasons. As the bivariate unidirectional causal methods, the CCF and TE methods miss the link I → S because the opposite direction S → I has a more significant causal effect. It should be noted that since the CCM method successfully identified all of the predefined causal links, many new false connections, i.e., Q → P , Q → S , S → P , I → P , were generated with the control of confounding and mediating variables. One reasonable explanation is that the synthetic system is strongly coupled, such that the CCM method easily mistakes unidirectional causal relationships for bidirectional relationships. In Model S2, PCMCI + perfectly restores all predefined causal links, while the other three methods (CCF, CCM, and TE) incorrectly identified many indirect causal links due to the transitivity of causal relationships. Similar to the results of Model 1, the CCF and TE methods miss the causal link I → S , and CCM also mistakes all unidirectional links for bidirectional ones. Table 3 further lists the evaluation results of the causal tests in the synthetic case. The PCMCI + method shows the best overall performance in both model structures, with the TPR, AR, and FPR values of 0.97(1.00), 0.85(0.96), and 0.23(0.05) in Model S1(S2). The CCF and TE methods show relatively lower TPR and AR values, and relatively higher FPR values. In contrast, the CCM method shows the lowest AR (lower than 0.5) and the highest TPR and FPR (higher than 0.9). Besides, in situations where direct and indirect causality are not distinguished, that is, the indirect links identified by the methods are not considered as the wrong links (the values in parentheses of Table 3 ), the performances of CCF, CCM, and TE methods improve significantly, especially in Model S2, indicating that the three methods are strongly affected by mediating variables and the identified causal links are likely to contain indirect effects, which might hinder our understanding of the real mechanism in hydrological systems. It should be noted that in such situations, the TPR values are not changed because the true positives (TP, i.e., the causal links generated by the hydrological model are correctly identified) remain the same. In addition, due to the complexity of the hydrological model, the overall performance of all causal methods in Model S2 is better than that in Model S1 except for the CCM method, which shows extremely high TPR and FPR simultaneously owing to the mistakes for bidirectional links. Therefore, in the following real case studies, the identified causal links were revised to unidirectional links, that is, the direction with a significant maximum causal effect was regarded as the only causal direction, thus to avoid misidentification. Table 3 Evaluation results of the causal tests in the synthetic case. CCF CCM TE PCMCI+ TPR 0.60 1.00 0.60 0.97 Model S1 AR 0.58 (0.60) 0.42 (0.50) 0.66 (0.69) 0.85 (0.87) FPR 0.43 (0.40) 1.00 (1.00) 0.30 (0.22) 0.23 (0.22) TPR 0.75 1.00 0.75 1.00 Model S2 AR 0.67 (0.89) 0.36 (0.49) 0.67 (0.89) 0.96 (1.00) FPR 0.38 (0.00) 0.95 (0.93) 0.38 (0.00) 0.05 (0.00) Note: The values in parentheses represent situations where direct and indirect causality are not distinguished, i.e., the indirect links identified by the methods are not considered as the wrong links. The TPR values are not changed in such situations. 3.3 Performance in sensitivity tests In this section, the sensitivity of each method in sample length, noise, and missing data was analyzed. Figure 4 presents the impact of sample length with the sizes of 100, 300, 500, 1000, and 2000. 100 datasets were generated to ensure the robustness of the results. As shown in Figs. 4 a and b, in both Model S1 and S2, for CCF and CCM methods, the TPR is not sensitive to variations in sample length; while for TE and PCMCI + methods, the TPR increases with increasing sample length. For the CCF method, which is based on linear lag correlation, 100 samples are enough to detect all the causal links. This number is also applicable to the CCM method, which is based on the deterministic dynamical system theory. In contrast, TE is based on the probability density framework, which needs to estimate the probability density function from the histogram of the frequency distribution, as well as to implement statistical hypothesis testing for conditional independence, thus requiring sufficient sample length. Similarly, the PCMCI + method also requires sufficient length samples in conditional independence tests to keep iterating and removing initialized redundant connections. To achieve relative stability, the TE and PCMCI + methods need at least 1000 and 500 samples for Model S1 and S2, respectively. This difference results from the complexity of the model structure. As for the Accuracy Rate (AR), the values in each causal method fluctuate or even decrease with the increasing sample length, especially for TE and CCM methods (Figs. 4 c and d). This is attributed to the limitations of causal methods in dealing with the impacts of confounding and mediating variables, which become more striking with the increase of simple length and lead to a significant increase in false positive rate. Additionally, the dashed lines in Fig. 4 represent the situations where direct and indirect causality are not distinguished. The difference in trends between the solid and dashed lines can be further analyzed to determine whether variations in sample length affect the identification of indirect causal links. The CCF, CCM, and PCMCI + methods show a similar trend between the solid and dashed lines in both model structures, while for the TE method, the variation rate is slightly inconsistent or even the reverse, due to the increasing detection rate of indirect links P → I in Model S1, and T → S , T → Q , M → Q in Model S2, with the increase of sample length. Besides, in the situation where direct and indirect causality is not distinguished, the improvement in the performance of causal methods in Model S2 is more pronounced than that in Model S1, which is consistent with the setup of model structure, i.e., the former is mainly controlled by indirect causality. Figure 5 presents the impact of noise with the SNR(dB) of 2, 3, 5, 10, 20, 30, and 40. The sample length was fixed as 2000, and 100 datasets were generated to ensure the robustness of the results. As shown in Fig. 5 a, in model S1, the noise has little impact on the TPR values for the CCF and CCM methods with the increase of noise level (i.e., the decrease of SNR). In contrast, for the TE method, the TPR increases with the increasing noise level, especially when the SNR (dB) is below 20. This is attributed to the assumption of the nonlinear stochastic system for the TE method, and appropriate noise helps to identify causal relationships. For PCMCI+, the TPR decreases with the decreasing SNR(dB) and shows a slight increase when the SNR (dB) is lower than 10. However, in Model S2, the noise has little impact on the TPR values for all causal methods (Fig. 5 b), owing to the insensitivity of this model structure to input noise. As for the Accuracy Rate (AR), the values in CCF and TE are relatively stable with the variations of the noise level in both model structures (Figs. 5 c and d), while for the PCMCI + method, the AR decreases with the decreasing SNR in Model S1, and for the CCM method, the AR increases as the SNR(dB) drops below 10 in both model structures. Additionally, all causal methods show a similar trend between the solid and dashed lines in both model structures, indicating that the noise does not have significant impacts on the identification of indirect causal links. Figure 6 presents the impact of missing data on the performance of each causal method with the sparsity rate increasing from 10–60%. Two scenarios of missing data, namely synchronous and asynchronous sparsity, which may occur in real hydrologic data due to equipment errors and defective storage technologies, were constructed with 100 tests. To ensure the computability of all causal methods, the missing values were complemented by the arithmetic mean of observed values of the same variable. As shown in Figs. 6 a and b, with the increase of synchronous missing rates, the TPR values remain relatively stable for the CCF method and begin to decline at a critical missing value for the other three methods, especially for CCM and PCMCI+. In comparison with the synchronous sparsity, the TPR values remain relatively stable for the PCMCI method in the scenario of asynchronous sparsity, while for the CCM method, the TPR values decline significantly over the 30% point of missing rate (Figs. 6 e and f). One reasonable explanation is that the CCM method requires state space reconstruction by continuous time series data, while the missing data can lose numerous effective information for the state space, especially in the asynchronous scenario, thus disturbing the detection of causal relationships. As for the Accuracy Rate (AR), with the increasing missing rate, the values in the CCF method remain relatively stable, while in the other three methods (CCM, TE, and PCMCI+), the values show fluctuation (Figs. 6 c, d, f and g). It should be noted that for the CCM method, the AR values increase substantially, especially over the critical missing value of 30%, which seems implausible. One reasonable explanation is that as the missing rate increases, effective information of variables in their state space is gradually lost, and causally linked variables in the system may no longer maintain the information signature of each other; therefore the efficiency of cross-mapping decreases. The detection rate of both correct and incorrect causal links declines simultaneously, and the latter is more pronounced. In addition, the CCF and PCMCI + methods show a similar trend between the solid and dashed lines, while for the TE method, the variation rate is slightly inconsistent or even the reverse over the 50% point of missing rate due to the decreasing detection rate of indirect links P → I in Model S1, and T → S and T → Q , M → Q in Model S2, with the increasing missing rate. For the CCM method, the variation rate is extremely inconsistent over the 30% point of the missing rate due to the declining detection rate of all indirect links ( P → I, I → Q, T → S , T → Q , M → Q ) simultaneously. At the end of this section, the results of the sensitivity tests can be summarized as follows: In the impact of sample length, the CCF and CCM methods show relatively stable performance, while the TE and PCMCI + methods require at least 500 samples to achieve relative stability; In the impact of noise, the performance of causal models is affected by the model structure, and the TE and PCMCI + methods perform more unstably for Model S1; In the impact of missing data, all causal methods show relatively stable performance within 30% of the missing rate, and the performance of the CCM method deteriorates rapidly over 30%. 4. Application to real case study In this section, two real study cases are presented to further evaluate the applicability of different causal inference methods in complex hydrological systems. 4.1 Application to real case 1 4.1.1 Study area and data The Shale Hills Catchment (SHC; Fig. 7 (a)), located in central Pennsylvania, USA, is a V-shaped small (0.08 km 2 ) forested headwater catchment with comparatively steep slopes and narrow ridges (Guo et al., 2014 ). The surface elevation ranges from 256 to 310 m, with relatively homogenous land cover and lithology, and relatively heterogeneous soil thickness and organic matter content (Lin, 2006 ). The catchment belongs to a temperate continental climate region and the average air temperature varies from − 3°C (January) to 22°C (July). The annual average precipitation is around 980 mm, and the monthly distribution is relatively uniform with a small maximum in summer when the rainfall is usually characterized with high intensity and short duration (Jiang et al., 2023 ). The first-order stream of the catchment converges to the Juniata River, and is usually dry during summer months (Liu et al., 2020 ). Besides, given that the snowfall mainly occurs from November to April, this period was set as the snow-cover period in this study, while the other period over the water year (May to October) was divided as the snow-free period, thus to explore different causal mechanisms during wet and dry conditions. The hydrological series, including discharge ( Q ), precipitation ( P ), groundwater level ( GW ), soil moisture ( SM ), and snow depth ( SD ), were obtained from the Critical Zone Observatory Data Site ( http://www.czo.psu.edu/data_time_series.html ). The dataset is a reanalysis result from the Flux-PIHM model (Shi et al., 2013 ), and has been subject to strict quality control. The temporal resolution is hourly with the period of 2009/11/1 ~ 2010/10/31. The min-max normalized values of the time series in snow-free (May to October) and snow-cover (November to April) periods are shown in Figs. 7 (b) and (c), respectively. 4.1.2 Results Figure 8 presents the causal structures identified by different causal inference methods for the Shale Hills Catchment (SHC) during the snow-free period. Generally, all methods can identify the main causal relationships P → SM , P → GW and P → Q , which represent the basic hydrological processes from precipitation to soil water, groundwater and streamflow, respectively, but show differences in other relationships. For the link between GW and Q , The CCF method presents a forward direction, while the CCM and TE methods present a backward direction. Only the PCMCI + method shows a bidirectional relationship, which contributes to the understanding of potential interaction processes between groundwater and streamflow. For the link between SM and GW , both CCF and CCM methods present a lagged forward link, while the PCMCI + method shows a contemporary backward link. Besides, the CCF and CCM methods present a forward link between SM and Q , which is omitted by the other two methods. In terms of causal strength, all methods support a relatively weak association for P → Q and P → GW , while the CCF and CCM methods tend to favor a stronger causal interaction between GW and Q , and the TE and PCMCI + methods tend to favor a stronger link P → SM. Figure 9 presents the causal structures identified by different causal inference methods for the SHC during the snow-cover period. The snow depth (SD) variable is added to this hydrological system. Similar to the results during the snow-free period, all methods can detect the main causal relationships P → Q and P → GW , and show differences in other relationships between GW and Q , SM and Q , and SM and GW . However, limitations of causal methods are exposed to this more complex hydrological system. The unreasonable causal links GW → SD and Q → SD identified by the CCF and CCM methods might be attributed to their shortcomings in dealing with confounding variables. The causal link SD → SM identified by CCF, CCM, and PCMCI + methods is consistent with the snowmelt process. However, distinguished from the negative causal strength in the PCMCI + method, the CCF method shows an incomprehensible positive causal strength, which might be attributed to the influence of the confounding variable P (simultaneously affecting SM and SD ). In comparison with the snow-free period, all causal methods in the snow-cover period support a weaker association between P and SM and a stronger association between P and GW except for the TE method (due to the lack of this causal link) during the snow-cover period. 4.2 Application to real case 2 4.2.1 Study area and data The Chuosijia River Basin (CRB; Fig. 10 (a)) is situated on the southeastern edge of the Tibetan Plateau, China, with a drainage area of 14813 km 2 and elevation ranging from 2440 to 5403m (Yang et al., 2023 ). Impacted by the westerly circulation and the southwest monsoon, the monsoon climate here is remarkable with distinct dry and wet seasons, and the same period of rain and heat. The multi-annual average temperature and precipitation are 8.6 ℃ and 740mm, respectively. The intra-annual distribution of precipitation is uneven. More than 80% of the precipitation is concentrated from June to October and usually in the form of heavy rainfall due to the warm moisture from the southwestern Indian Ocean (Chen and Alexander, 2022 ). The multi-year average runoff depth is 392mm, and runoff is primarily formed by precipitation, followed by snowmelt and groundwater. The first-order stream of the basin finally converges to Minjiang River, a main tributary of the upper Yangtze River (Liang et al., 2023 ). Besides, the CRB has little human interference and can be deemed as a natural basin, which is suitable for investigating causal interactions in natural hydrological systems. Considering that the snowfall mainly occurs from November to May, this period was set as the snow-cover period in this study and the other period over the water year (June to October) was divided as the snow-free period, thus to explore different causal mechanisms during dry and wet conditions. The hydrological series includes precipitation ( P ), soil moisture ( SM ), evaporation ( E ), snow water equivalent ( SWE ), and discharge ( Q ). The precipitation data was obtained from CN05.1(Wu and Gao, 2013 ), a gridded dataset (0.25°×0.25°; https://ccrc.iap.ac.cn/ ) based on more than 2400 observed meteorological stations. The soil moisture data (0-100 cm) was gained from SMCI1.0 (Li et al., 2022 ), a long-term high-resolution soil moisture dataset (1km×1km; https://data.tpdc.ac.cn/home/ ) based on the measurements of 1789 observed stations across China. The evaporation and snow water equivalent data were obtained from GLEAM (0.25°×0.25°; https://www.gleam.eu/ ) and ERA5-land (0.1°×0.1°; https://cds.climate.copernicus.eu/#!/home/ ), respectively, which have been widely applied in the Tibetan Plateau recently (Chen et al., 2022 ; Li et al., 2019 ; Yang et al., 2020 ). The discharge data was gained from the Hydrological Yearbook of the Yangtze River Basin. All gridded data was aggregated to the whole basin scale, with the daily temporal resolution and the whole period from 2008 to 2018. The min-max normalized values of the time series in snow-free (June to October) and snow-cover (November to May) periods during water year 2015–2016 are shown in Figs. 10 (b) and (c), respectively. 4.2.2 Results Figure 10 presents the causal structures identified by different causal inference methods for the Chuosijia River Basin (CRB) during the snow-free period. In general, all methods can identify the main causal relationships P → SM and P → Q , which denote the basic hydrological processes from precipitation to soil water and streamflow, respectively, but show differences in other relationships. For the link between SM and Q , the CCF and CCM methods present a lagged forward link, while the PCMCI + method shows a contemporary backward link. For the link between E and P , all methods show a forward direction except for the CCM method. For the links between E and SM , and Q and E , the CCF and PCMCI + methods present forward links, while the CCM method shows backward links, providing a complementary understanding of the interactions among daily evaporation, soil moisture, and streamflow. With respect to the causal strength, all methods support a relatively strong association for P → Q. Figure 11 shows the causal structures identified by different causal inference methods for the CRB during the snow-cover period. The snow water equivalent ( SWE ) variable is added to this hydrological system. Generally, all methods can identify the causal relationships E → P and P → Q , which reflect the basic processes of the hydrologic cycle, i.e., from evaporation to precipitation and then to streamflow, but show large differences in other relationships. For the link between P and SM , the CCF and PCMCI + methods present a forward link, while the CCM and TE methods show a backward link. Besides, the CCF and CCM methods identify the causal links E → Q and E → SM , which are omitted by the other two methods. Besides, except for the TE method, all methods identify the causal links SM → Q and E → SWE , which reflect subsurface flow and snowmelt processes, respectively. Only the CCM method detects the snow cover process ( P → SWE ). Nevertheless, many unreasonable causal relationships emerge in this five-variable system, for instance, the SWE → P detected by the CCF and TE methods, and the Q → SWE detected by the CCM method. Compared with the snow-free period, all causal methods in the snow-cover period support a stronger association between E and P and a weaker association between P and Q in the snow-cover period, and the CCF and CCM methods support a stronger association between SM and Q . The results indicates the increasing contribution of baseflow to streamflow during the snow-free period, as well as the diminishing effect of precipitation on streamflow. Discussion 5.1 Comparison of four causal inference methods Identifying causal associations in hydrological systems is challenging. Applying different causal methods to artificially generated and real-world hydrological data has been confirmed to yield some incomprehensible or even contradictory results, which leads to some reasonable doubts about the reliability of the currently popular methods, and highlights the importance of maintaining critical attention to the limitations of individual methods. The CCF method, which is rooted in traditional linear lag correlation theory, shows numerous significant connections whether in synthetic or real cases. Since reducing the significance level may control the number of false links (Rinderera et al., 2018 ), yet, significant associations can still be detected as the p -value decreases to 0.001. This might be attributed to relatively strong coupling and synchronization of the hydrological systems. Besides, the essence of the CCF method lies in statistical dependencies, without excluding indirect and confounding factors (Dean and Dunsmuir, 2016 ), which leads to many spurious causal associations, such as Q → S (Fig. 3 a), M → I (Fig. 3 e) in the synthetic case, and GW → SD (Fig. 9 a), SWE → P (Fig. 12 a) in the real cases 1 and 2, respectively. Nevertheless, the CCF method is simple and efficient, and presents a stable performance in the sensitivity tests, with little variations induced by the sample length, noise, and missing values. Therefore, the CCF method can be used as a preliminary means of inferring causal associations in hydrological systems, and a benchmark for the comparison of other causal inference methods. For the CCM method, results in the synthetic case show extremely high false positive rates (almost 1) in both model structures, even if all preset causal relationships were identified simultaneously (Figs. 3 b, f). The irrational bi-directional causality was also reported by Ombadi et al. ( 2020 ), which highlights the difficulty for the CCM method to address strong coupling in hydrological systems. In real cases, the causal directions were revised to that with the significant maximum causal effect, while some unreasonable causal links were still identified, such as GW → SD , Q → SD (Fig. 9 b), and Q → SWE (Fig. 12 b) in the real case 1 and 2, respectively, which further confirms the limitations of CCM in dealing with confounding factors (Delforge et al., 2022 ). Nevertheless, the CCM method still helps to understand complex hydrological processes, which requires us to return to the basic principles and assumptions of this method, that is, the nonlinear causality identified by CCM would only make sense under the framework of nonlinear dynamic system (NDS) with its concept fully understood (Zhao et al., 2023a ). The CCM method might contribute to exploring hidden hydrological causal mechanisms from the perspective of NDS, such as the interactions between groundwater and streamflow (Bonotto et al., 2022 ), and within cryospheric hydrological cycle (Zhao et al., 2023b ), which may not be directly observed from the geophysical level, and thus in turn helps to develop models based on physical foundations. For the TE method, the results exhibit better performance in synthetic large sample tests, compared with the CCF and CCM methods, especially in cases where direct and indirect causality are not distinguished. In real cases, the identified significant causal relationships are relatively few, but most of them can be explained by basic hydrological laws, such as the links P → SM (Figs. 8 c, 9 c) and P → Q (Figs. 11 c, 12 c). Based on information theory, the TE method is not constrained by linear assumptions and partly overcomes the autocorrelation of hydrological time series using conditional mutual information, thus can be well applied in complex hydrological systems even with the threshold effect (Moges et al., 2022 ; Tennant et al., 2020 ). However, similar to the CCF and CCM methods, it fails to fundamentally avoid the interference of confounders and indirect causality, and some unreasonable connections still need to be cautiously explained. Besides, the sensitivity analysis indicates the necessity of sufficient sample length (i.e., at least 500–1000) for TE analysis, due to the challenge of estimating three-dimensional probability density, which has also been emphasized in other research (Kathpalia et al., 2022 ; Ruddell and Kumar, 2009 ) The PCMCI + method, which is based on causal network learning, shows optimal performance in synthetic large-sample tests, and satisfactory interpretability in real cases. This method infers causality from the network graph of multivariate time series, utilizing iterative conditional independence tests to largely overcome the interference of autocorrelation and confounding factors, and distinguish direct and indirect causality in hydrological systems (Fig. 3 d, h). Notwithstanding, contemporary causal associations in real cases, partly due to the coarse time resolutions (Runge et al., 2019a ), need to be explained with caution. In terms of the algorithms for the conditional independence test, partial correlation (ParCorr) was selected in this study, which has also been widely adopted in other research (Karmouche et al., 2023 ; Nowack et al., 2020 ). Since other testing algorithms, such as conditional mutual information (CMI), seems suitable for discovering nonlinear associations in hydrological systems, yet it is fairly time-consuming (Runge et al., 2019b ), and recent research showed a low recall rate in synthetic datasets, and an unstable performance in real cases (Delforge et al., 2022 ), which highlights the importance of choosing appropriate testing algorithms. Moreover, the graph-based method requires the assumption of causal sufficiency, i.e. the absence of hidden variables. In general, this assumption cannot be directly tested, and researchers can only consider known factors, as many as possible, in the causal network analysis through their prior expert knowledge, which may add potential uncertainty to the results. Yet, some recent work indicates that the assumption could be appropriately relaxed, especially in the real world (Gerhardus, 2020 ; Runge, 2018 ). In summary, we recommend a flexible strategy for selecting causal inference methods based on the purpose of research. For exploratory work, such as investigating possible causal connections, the CCF and CCM methods are suggested, even if some irrational links are inevitably produced. The former detects all possible linear associations in hydrological systems, while the latter can reveal the hidden causal interactions in nonlinear dynamic systems. For confirmatory work, such as selecting significant predictors, the TE method can be selected with relatively low false detections. The PCMCI + method, with optimal integrative performance and remarkable interpretability, can be applied to both works mentioned above, and is suggested to deeply understand the interaction mechanisms among multivariates in complex hydrological systems. Nevertheless, each method improves, more or less, our understanding of hydrological processes, which calls into our multi-perspective comprehension of the findings obtained by these methods under the context of their particular assumptions and constraints. Specifically, the causality identified by the CCF method can be understood as the linear connection, while identified by the CCM, TE, and PCMCI + methods need to be understood from nonlinear dynamics, information theory, and causal network perspectives, respectively. Therefore, it is important to combine the priori expert knowledge with the assumptions and limitations of each method, to comprehensively explain the reasonable or implausible causal associations. Moreover, the application of different causal methods may yield conflicting results, which highlights the importance to state the assumptions when drawing conclusions of causal issues. To this end, hydrological researchers are encouraged to be more explicit in elaborating assumptions that enable more robust conclusions, and in interpreting and evaluating conclusions under alternative sets of assumptions (Runge et al., 2023 ). 5.2 Direct and indirect causality As the central focus of hydrology, the water cycle contains numerous components, including rainfall, evaporation, soil moisture, runoff, groundwater, snow, etc., which are connected in direct (e.g., rainfall-soil moisture) or indirect (e.g. groundwater-evaporation) way. Causal inference contributes to revealing the functional connectivity of these components in hydrological systems (Delforge et al., 2022 ; Rinderera et al., 2018 ). However, the causal methods may erroneously identify indirect causal influences as direct ones due to the transitivity of the causal relationship (Park et al., 2023 ). On the other hand, if two components are not connected directly, a causal influence may also exist, but it must be indirect. It is of great importance to remove the effects of indirect causality and thus to determine the direct causal relationships, as the latter serves as the foundation for modeling, prediction, and control of the system (Leng et al., 2020 ). In comparison with other algorithms, the PCMCI + method can effectively distinguish between direct and indirect causal links (Fig. 3 ; Table 3 ), with explicit indirect impact mechanisms from the graph of multivariate time series. Nevertheless, discriminating between direct and indirect causality for the real hydrological systems is still challenging, because both relationships are often intertwined. Taking the P - SM-Q process for example, the precipitation falling on land can first supplement soil moisture through infiltration and then move to the outlet of the basin through subsurface processes (i.e., the indirect process), or directly transfer through surface processes such as saturation overland flow or infiltration excess runoff (i.e., the direct process) (Kidron, 2021 ). How to quantify such complex causal interactions among three or even more variables in hydrological systems remains to be further considered (Goodwell and Bassiouni, 2022 ; Weijs et al., 2018 ). 5.3 Progress, limitations, and future perspectives In comparison with the previous research (Delforge et al., 2022 ; Ombadi et al., 2020 ), our study focuses on hydrological systems, introduces the latest causal network learning algorithm PCMCI+, and serves the CCF method as the link between causality and correlation, and as the benchmark for the comparison of other methods. Moreover, in the synthetic cases, we expanded the structures of the conceptual hydrological models, enriched the sensitivity tests, and discussed the direct and indirect causality in hydrological systems for the first time. Two real cases deepen our understanding of different causal methods at different spatiotemporal scales, respectively. Our research in both cases systematically investigated the question of when and how to apply these methods, and how to interpret their results, thus contributing to their further popularity in the hydrology community. Nonetheless, this study concentrates more on the methodological issues, the inner mechanisms of hydrological processes across spatio-temporal scales in the study area are beyond the scope of this study, which need further research (Liu et al., 2020 ; Wen et al., 2021 ; Hao et al., 2022 ). Besides, under different causal frameworks, testing causal mechanisms either unduly relies on assumptions or lacks theoretical examination (Su et al., 2023 ). Thus, the discovery and test of causal mechanisms in complicated hydrological systems with unknown causal structures is still challenging. Moreover, the hydrological systems may involve certain latent variables, such as unknown/hidden or unobservable variables, which can introduce some confounding factors and make the detected links spurious. Fortunately, the LPCMCI algorithm, proposed by Runge ( 2020 ) recently, may overcome this limitation to some extent and can be further applied in hydrology. Additionally, some new methods, such as the Partial cross mapping (PCM; Leng et al., 2020 ) and the method proposed by Park et al. ( 2023 ), which address the issue of indirect causality based on nonlinear dynamics theory and monotonic ODE model, respectively, can be tested and evaluated in the future. Recently, the questions–assumptions–data (QAD) template, proposed by Runge et al. ( 2023 ), helps to guide researchers on how to phrase and tackle their issues in the framework of causal inference. Yet, some typical challenges, such as the hidden confounding, non-stationarity, contemporaneous causality, and preprocessing of variables, need extra attention. Fortunately, the era of big data provides new opportunities for research in this field, as it could comprehensively characterize the structural information between variables, verify the spurious causal links in the causal network structure, and infer latent variable structures that are difficult to observe (Li et al., 2023). To this end, causal inference helps to build the bridge between data-driven machine learning and prior expert knowledge, thus promoting the understanding of complex hydrological processes and the robust causal prediction. Conclusion In this research, the performances of four popular causal inference methods (CCF, CCM, TE, and PCMCI+) were systematically evaluated in hydrological systems using both synthetic and real-world time series. For the synthetic cases, the PCMCI + method shows the best performance in large sample tests, while the other three methods present relatively poor performances due to the interference of confounding and mediating factors. In sensitivity tests, the CCF method is less affected by sample length, noise, and missing values, while the CCM method is significantly impacted by the missing rate, especially over the critical value of 30%. Besides, the TE and PCMCI + methods need at least 500 samples to achieve relative stability, and are vulnerable to the interference of noise. For the real cases, the PCMCI + method shows the best interpretability, while the other three methods produce many inexplicable causal links. Additionally, some contradictory results among different methods emerge, owing to their different assumptions and limitations. In summary, the PCMCI + method serves as a favorable choice for conducting causal inference within hydrological systems, while some strong assumptions, such as the causal sufficiency, need to be considered with caution. Nonetheless, each method improves, more or less, our understanding of hydrological processes, which requires our multi-perspective comprehension of the findings obtained by these methods under the context of their particular assumptions and constraints. A comprehensive application of diverse methods based on the specific issue is encouraged for the robustness of conclusions, with the assumptions stated explicitly in advance. Promisingly, the causal inference methods provide a complementary data-driven avenue to unlock the inner mechanisms of complex hydrological systems, and have broad application prospects in the hydrology community. Declarations Acknowledgements The research is supported by the National Nature Science Foundation of China (No. 52379017). CRediT authorship contribution statement Hanxu Liang: Methodology, Investigation, Formal analysis, Visualization, Writing-original draft. Wensheng Wang: Supervision, Funding acquisition, Writing - Review & Editing. Bin Chen: Data curation. Li Guo: Writing - Review & Editing. Hu Liu: Writing - Review & Editing. Siyi Yu: Validation. Dan Zhang: Conceptualization, Validation. Declaration of Competing Interest All authors agreed to the published version of the manuscript and declare no conflicts of interests. References Altman, N., and M. Krzywinski (2015), Association, correlation and causation, Nat. Methods 12 (10), 899-900. https://doi.org/10.1038/nmeth.3587. Angelov, P. P., E. A. Soares, R. C. Jiang, N. I. Arnold, and P. M. Atkinson (2021), Explainable artificial intelligence: an analytical review, Wires. Data. Min. Knowl. 11 (5). https://doi.org/10.1002/widm.1424. Apaydin, H., and M. Sibtain (2021), A multivariate streamflow forecasting model by integrating improved complete ensemble empirical mode decomposition with additive noise, sample entropy, Gini index and sequence-to-sequence approaches, J. Hydrol. 603 . https://doi.org/10.1016/j.jhydrol.2021.126831. Bennett, A., B. Nijssen, G. X. Ou, M. Clark, and G. Nearing (2019), Quantifying Process Connectivity With Transfer Entropy in Hydrologic Models, Water Resour. Res. 55 (6), 4613-4629. https://doi.org/10.1029/2018wr024555. Bergström, S., and L. P. Graham (1998), On the scale problem in hydrological modelling, J. Hydrol. 211 (1-4), 253-265. https://doi.org/10.1016/s0022-1694(98)00248-0. Blöschl, G., M. F. P. Bierkens, A. Chambel, C. Cudennec, and G. Destouni (2019), Twenty-three unsolved problems in hydrology (UPH) - a community perspective, Hydrol. Sci. J. 64 (10), 1141-1158. https://doi.org/10.1080/02626667.2019.1620507. Bonotto, G., T. J. Peterson, K. Fowler, and A. W. Western (2022), Identifying Causal Interactions Between Groundwater and Streamflow Using Convergent Cross-Mapping, Water Resour. Res. 58 (8). https://doi.org/10.1029/2021wr030231. Budakoti, S., T. Chauhan, R. Murtugudde, S. Karmakar, and S. Ghosh (2021), Feedback From Vegetation to Interannual Variations of Indian Summer Monsoon Rainfall, Water Resour. Res. 57 (5). https://doi.org/10.1029/2020wr028750. Chen, J., F. P. Brissette, and R. Leconte (2010), A daily stochastic weather generator for preserving low-frequency of climate variability, J. Hydrol. 388 (3-4), 480-490. https://doi.org/10.1016/j.jhydrol.2010.05.032. Chen, F., W. T. Crow, P. J. Starks, and D. N. Moriasi (2011), Improving hydrologic predictions of a catchment model via assimilation of surface soil moisture, Adv. Water Resour. 34 (4), 526-536. https://doi.org/10.1016/j.advwatres.2011.01.011. Chen, R., M. X. Yang, X. J. Wang, G. N. Wan, and H. Y. Li (2022), Thermal regime variations of the uppermost soil layer in the central Tibetan Plateau, Catena 213 . https://doi.org/10.1016/j.catena.2022.106224. Chen, Y., and D. Alexander (2022), Integrated flood risk assessment of river basins: Application in the Dadu river basin, China, J. Hydrol. 613 . https://doi.org/10.1016/j.jhydrol.2022.128456. Cryer, J. D., and K. Chan (2008), Time series analysis with applications in R. New York, NY: Springer. Dean, R. T., and W. T. M. Dunsmuir (2016), Dangers and uses of cross-correlation in analyzing time series in perception, performance, movement, and neuroscience: The importance of constructing transfer function autoregressive models, Behav. Res. Methods 48 (2), 783-802. https://doi.org/10.3758/s13428-015-0611-2. Delforge, D., O. de Viron, M. Vanclooster, M. Van Camp, and A. Watlet (2022), Detecting hydrological connectivity using causal inference from time series: synthetic and real karstic case studies, Hydrol. Earth Syst. Sci. 26 (8), 2181-2199. https://doi.org/10.5194/hess-26-2181-2022. Faybishenko, B. (2017), Detecting dynamic causal inference in nonlinear two-phase fracture flow, Adv. Water Resour. 106 , 111-120. https://doi.org/https://doi.org/10.1016/j.advwatres.2017.02.011. Gao, Y. B., C. Merz, G. Lischeid, and M. Schneider (2018), A review on missing hydrological data processing, Environ. Earth. Sci. 77 (2). https://doi.org/10.1007/s12665-018-7228-6. Geng, S., F. Devries, and I. Supit (1986), A simple method for generating daily rainfall data, Agric. Meteorol. 36 (4), 363-376. https://doi.org/10.1016/0168-1923(86)90014-6. Gerhardus, A. a. R., Jakob (2020), High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems, volume 33, pages 12615–12625. Curran Associates, Inc. Gong, W., D. W. Yang, H. V. Gupta, and G. Nearing (2014), Estimating information entropy for hydrological data: One-dimensional case, Water Resour. Res. 50 (6), 5003-5018. https://doi.org/10.1002/2014wr015874. Good, S. P., D. Noone, and G. Bowen (2015), Hydrologic connectivity constrains partitioning of global terrestrial water fluxes, Science 349 (6244), 175-177. https://doi.org/10.1126/science.aaa5931. Goodwell, A. E., and M. Bassiouni (2022), Source Relationships and Model Structures Determine Information Flow Paths in Ecohydrologic Models, Water Resour. Res. 58 (9). https://doi.org/10.1029/2021wr031164. Granger, C. W. (1969), Investigating causal relations by econometric models and cross‐spectral methods, Econometrica: Journal of the Econometric Society, 37(3), 424–438 . https://doi.org/10.2307/1912791. Guo, L., J. Chen, and H. Lin (2014), Subsurface lateral preferential flow network revealed by time-lapse ground-penetrating radar in a hillslope, Water Resour. Res. 50 (12), 9127-9147. https://doi.org/10.1002/2013wr014603. Hao, Y., F. B. Sun, H. Wang, W. B. Liu, Y. J. Shen, Z. Li, and S. J. Hu (2022), Understanding climate-induced changes of snow hydrological processes in the Kaidu River Basin through the CemaNeige-GR6J model, Catena 212 . https://doi.org/10.1016/j.catena.2022.106082. Jiang, Y. J., Y. L. Zhang, B. H. Fan, J. H. Wen, H. Liu, C. R. Mello, J. F. Cui, C. Yuan, and L. Guo (2023), Preferential flow influences the temporal stability of soil moisture in a headwater catchment, Geoderma 437 . https://doi.org/10.1016/j.geoderma.2023.116590. Karmouche, S., E. Galytska, J. Runge, G. A. Meehl, A. S. Phillips, K. Weigel, and V. Eyring (2023), Regime-oriented causal model evaluation of Atlantic-Pacific teleconnections in CMIP6, Earth System Dynamics 14 (2), 309-344. https://doi.org/10.5194/esd-14-309-2023. Kathpalia, A., P. Manshour, and M. Palus (2022), Compression complexity with ordinal patterns for robust causal inference in irregularly sampled time series, Sci. Rep. 12 (1). https://doi.org/10.1038/s41598-022-18288-4. Khatibi, R., B. Sivakumar, M. A. Ghorbani, O. Kisi, K. Kocak, and D. F. Zadeh (2012), Investigating chaos in river stage and discharge time series, J. Hydrol. 414 , 108-117. https://doi.org/10.1016/j.jhydrol.2011.10.026. Kidron, G. J. (2021), Comparing overland flow processes between semiarid and humid regions: Does saturation overland flow take place in semiarid regions?, J. Hydrol. 593 . https://doi.org/10.1016/j.jhydrol.2020.125624. Kinney, J. B., and G. S. Atwal (2014), Equitability, mutual information, and the maximal information coefficient, Proc. Natl. Acad. Sci. U. S. A. 111 (9), 3354-3359. https://doi.org/10.1073/pnas.1309933111. Kleidon, A., and M. Renner (2013), Thermodynamic limits of hydrologic cycling within the Earth system: concepts, estimates and implications, Hydrol. Earth Syst. Sci. 17 (7), 2873-2892. https://doi.org/10.5194/hess-17-2873-2013. Konapala, G., S. C. Kao, and N. Addor (2020), Exploring Hydrologic Model Process Connectivity at the Continental Scale Through an Information Theory Approach, Water Resour. Res. 56 (10). https://doi.org/10.1029/2020wr027340. Leng, S. Y., H. F. Ma, J. Kurths, Y. C. Lai, W. Lin, K. Aihara, and L. N. Chen (2020), Partial cross mapping eliminates indirect causal influences, Nat. Commun. 11 (1). https://doi.org/10.1038/s41467-020-16238-0. Li, Q. L., G. S. Shi, W. Shangguan, V. Nourani, J. D. Li, L. Li, F. N. Huang, Y. Zhang, C. Y. Wang, D. G. Wang, J. X. Qiu, X. J. Lu, and Y. J. Dai (2022), A 1 km daily soil moisture dataset over China using in situ measurement and machine learning, Earth. Syst. Sci. Data 14 (12), 5267-5286. https://doi.org/10.5194/essd-14-5267-2022. Li, X., M. Feng, Y. H. Ran, Y. Su, F. Liu, C. L. Huang, H. F. Shen, Q. Xiao, J. B. Su, S. W. Yuan, and H. D. Guo (2023a), Big Data in Earth system science and progress towards a digital twin, Nat. Rev. Earth Env. https://doi.org/10.1038/s43017-023-00409-w. Li, Z. W., X. L. Xu, and K. L. Wang (2023b), Effects of distribution patterns of karst landscapes on runoff and sediment yield in karst watersheds, Catena 223 . https://doi.org/10.1016/j.catena.2023.106947. Li, X. Y., D. Long, Z. Y. Han, B. R. Scanlon, Z. L. Sun, P. F. Han, and A. Z. Hou (2019), Evapotranspiration Estimation for Tibetan Plateau Headwaters Using Conjoint Terrestrial and Atmospheric Water Balances and Multisource Remote Sensing, Water Resour. Res. 55 (11), 8608-8630. https://doi.org/10.1029/2019wr025196. Liang, H. X., D. Zhang, W. S. Wang, S. Y. Yu, and S. Nimai (2023), Evaluating future water security in the upper Yangtze River Basin under a changing environment, Sci. Total Environ. 889 . https://doi.org/10.1016/j.scitotenv.2023.164101. Lin, H. (2006), Temporal stability of soil moisture spatial pattern and subsurface preferential flow pathways in the shale hills catchment, Vadose Zone J. 5 (1), 317-340. https://doi.org/10.2136/vzj2005.0058. Liu, H., Y. Yu, W. Z. Zhao, L. Guo, J. T. Liu, and Q. Y. Yang (2020), Inferring Subsurface Preferential Flow Features From a Wavelet Analysis of Hydrological Signals in the Shale Hills Catchment, Water Resour. Res. 56 (11). https://doi.org/10.1029/2019wr026668. Luo, M., F. Meng, Y. Wang, C. Sa, Y. Duan, Y. Bao, and T. Liu (2023), Quantitative detection and attribution of soil moisture heterogeneity and variability in the Mongolian Plateau, J. Hydrol. 621 . https://doi.org/10.1016/j.jhydrol.2023.129673. Massei, N., J. P. Dupont, B. J. Mahler, B. Laignel, M. Fournier, D. Valdes, and S. Ogier (2006), Investigating transport properties and turbidity dynamics of a karst aquifer using correlation, spectral, and wavelet analyses, J. Hydrol. 329 (1), 244-257. https://doi.org/https://doi.org/10.1016/j.jhydrol.2006.02.021. Moges, E., B. L. Ruddell, L. Zhang, J. M. Driscoll, and L. G. Larsen (2022), Strength and Memory of Precipitation's Control Over Streamflow Across the Conterminous United States, Water Resour. Res. 58 (3). https://doi.org/10.1029/2021wr030186. Nowack, P., J. Runge, V. Eyring, and J. D. Haigh (2020), Causal networks for climate model evaluation and constrained projections, Nat. Commun. 11 (1). https://doi.org/10.1038/s41467-020-15195-y. Ombadi, M., P. Nguyen, S. Sorooshian, and K. l. Hsu (2020), Evaluation of Methods for Causal Discovery in Hydrometeorological Systems, Water Resour. Res. 56 (7). https://doi.org/10.1029/2020wr027251. Park, S. H., S. Ha, and J. K. Kim (2023), A general model-based causal inference method overcomes the curse of synchrony and indirect effect, Nat. Commun. 14 (1). https://doi.org/10.1038/s41467-023-39983-4. Patil, S., and M. Stieglitz (2014), Modelling daily streamflow at ungauged catchments: what information is necessary?, Hydrol. Process. 28 (3), 1159-1169. https://doi.org/10.1002/hyp.9660. Pearl, J. (2009), Causality: Models, Reasoning, and Inference, Cambridge University Press, Cambridge, 2 edn. https://doi.org/10.1017/CBO9780511803161. Peng, S. L., K. Mihara, X. L. Xu, K. Kuramochi, Y. Toma, and R. Hatano (2024), Modeling hydrological processes under Multi-Model projections of climate change in a cold region of Hokkaido, Japan, Catena 234 . https://doi.org/10.1016/j.catena.2023.107605. Reichenbach, H. (1956), The direction of time, University of California Press, Berkeley and Los Angeles . Reichstein, M., G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, and Prabhat (2019), Deep learning and process understanding for data-driven Earth system science, Nature 566 (7743), 195-204. https://doi.org/10.1038/s41586-019-0912-1. Richardson, C. W. (1981), Stochastic simulation of daily precipitation, temperature, and solar-radiation, Water Resour. Res. 17 (1), 182-190. https://doi.org/10.1029/WR017i001p00182. Rinderera, M., G. Ali, and L. G. Larsen (2018), Assessing structural, functional and effective hydrologic connectivity with brain neuroscience methods: State-of-the-art and research directions, Earth-Sci. Rev. 178 , 29-47. https://doi.org/10.1016/j.earscirev.2018.01.009. Ruddell, B. L., and P. Kumar (2009), Ecohydrologic process networks: 1. Identification, Water Resour. Res. 45 (3). https://doi.org/10.1029/2008wr007279. Rudin, C. (2019), Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead, Nat. Mach. Intell. 1 (5), 206-215. https://doi.org/10.1038/s42256-019-0048-x. Runge, J. (2018), Causal network reconstruction from time series: From theoretical assumptions to practical estimation, Chaos 28 (7). https://doi.org/10.1063/1.5025050. Runge, J. (2020), Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Proc. 36th Conf. Uncertainty in Artificial Intelligence (UAI) Vol. 124 of Proc. Machine Learning Research (eds Peters, J. & Sontag, D.) 1388–1397. . Runge, J., S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Munoz-Mari, E. H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Scholkopf, P. Spirtes, G. Sugihara, J. Sun, K. Zhang, and J. Zscheischler (2019a), Inferring causation from time series in Earth system sciences, Nat. Commun. 10 . https://doi.org/10.1038/s41467-019-10105-3. Runge, J., A. Gerhardus, G. Varando, V. Eyring, and G. Camps-Valls (2023), Causal inference for time series, Nat. Rev. Earth Env. 4 (7), 487-505. https://doi.org/10.1038/s43017-023-00431-y. Runge, J., P. Nowack, M. Kretschmer, S. Flaxman, and D. Sejdinovic (2019b), Detecting and quantifying causal associations in large nonlinear time series datasets, Sci. Adv. 5 (11). https://doi.org/10.1126/sciadv.aau4996. Sang, Y. F., V. P. Singh, J. Wen, and C. M. Liu (2015), Gradation of complexity and predictability of hydrological processes, J. Geophys. Res. Atmos. 120 (11), 5334-5343. https://doi.org/10.1002/2014jd022844. Schreiber, T. (2000), Measuring information transfer, Phys. Rev. Lett. 85 (2), 461-464. https://doi.org/10.1103/PhysRevLett.85.461. Shannon, C. E. (1948), A mathematical theory of communication, Bell System Technical Journal, 27(3), 379–423 . https://doi.org/10.1002/j.1538-7305.1948.tb01338.x. Shi, Y., K. J. Davis, C. J. Duffy, and X. Yu (2013), Development of a Coupled Land Surface Hydrologic Model and Evaluation at a Critical Zone Observatory, J. Hydrometeorol. 14 (5), 1401-1420. https://doi.org/10.1175/jhm-d-12-0145.1. Shojaie, A., and E. B. Fox (2022), Granger Causality: A Review and Recent Advances, Annu. Rev. Stat. Appl. 9 , 289-319. https://doi.org/10.1146/annurev-statistics-040120-010930. Sivakumar, B. (2004), Chaos theory in geophysics: past, present and future, Chaos Solitons & Fractals 19 (2), 441-462. https://doi.org/10.1016/s0960-0779(03)00055-9. Student (1908), The Probable Error of a Mean, Biometrika 6 (1), 1-25. https://doi.org/10.2307/2331554. Su, J. B., D. X. Chen, D. H. Zheng, Y. Su, and X. Li (2023), The insight of why: Causal inference in Earth system science, Science China-Earth Sciences . https://doi.org/10.1007/s11430-023-1148-7. Sugihara, G., E. R. Deyle, and H. Ye (2017), Misconceptions about causation with synchronyand seasonal drivers reply, Proc. Natl. Acad. Sci. U. S. A. 114 (12), E2272-E2274. https://doi.org/10.1073/pnas.1700998114. Sugihara, G., R. May, H. Ye, C. H. Hsieh, E. Deyle, M. Fogarty, and S. Munch (2012), Detecting Causality in Complex Ecosystems, Science 338 (6106), 496-500. https://doi.org/10.1126/science.1227079. Takens, F. (1981), Detecting Strange Attractors in Turbulence, In Dynamical systems and turbulence, Warwick 1980, (pp. 366–381). Berlin, Heidelberg: Springer . https://doi.org/10.1007/BFb0091924. Tennant, C., L. Larsen, D. Bellugi, E. Moges, L. Zhang, and H. X. Ma (2020), The Utility of Information Flow in Formulating Discharge Forecast Models: A Case Study From an Arid Snow-Dominated Catchment, Water Resour. Res. 56 (8). https://doi.org/10.1029/2019wr024908. Wang, Q. J., C. F. Yue, X. Q. Li, P. Liao, and X. Y. Li (2023), Enhancing robustness of monthly streamflow forecasting model using embedded-feature selection algorithm based on improved gray wolf optimizer, J. Hydrol. 617 . https://doi.org/10.1016/j.jhydrol.2022.128995. Wang, Y., J. Yang, Y. Chen, P. De Maeyer, Z. Li, and W. Duan (2018), Detecting the Causal Effect of Soil Moisture on Precipitation Using Convergent Cross Mapping, Sci. Rep. 8 . https://doi.org/10.1038/s41598-018-30669-2. Weijs, S. V., H. Foroozand, and A. Kumar (2018), Dependency and Redundancy: How Information Theory Untangles Three Variable Interactions in Environmental Data, Water Resour. Res. 54 (10), 7143-7148. https://doi.org/10.1029/2018wr022649. Wen, H., S. L. Brantley, K. J. Davis, J. M. Duncan, and L. Li (2021), The Limits of Homogenization: What Hydrological Dynamics can a Simple Model Represent at the Catchment Scale?, Water Resour. Res. 57 (6). https://doi.org/10.1029/2020wr029528. Wu, J., and X. J. Gao (2013), A gridded daily observation dataset over China region and comparison with the other datasets, Chinese Journal of Geophysics-Chinese Edition 56 (4), 1102-1111. https://doi.org/10.6038/cjg20130406. Xu, Z. P., X. L. Man, L. L. Duan, and T. J. Cai (2022), Improved subsurface soil moisture prediction from surface soil moisture through the integration of the (de)coupling effect, J. Hydrol. 608 , 12. https://doi.org/10.1016/j.jhydrol.2022.127634. Yang, W. J., Y. B. Wang, X. Liu, H. P. Zhao, R. Shao, and G. X. Wang (2020), Evaluation of the rescaled complementary principle in the estimation of evaporation on the Tibetan Plateau, Sci. Total Environ. 699 . https://doi.org/10.1016/j.scitotenv.2019.134367. Yang, Y., S. J. Chen, Y. R. Zhou, G. W. Ma, W. B. Huang, and Y. M. Zhu (2023), Method for quantitatively assessing the impact of an inter-basin water transfer project on ecological environment-power generation in a water supply region, J. Hydrol. 618 . https://doi.org/10.1016/j.jhydrol.2023.129250. Ye, H., E. R. Deyle, L. J. Gilarranz, and G. Sugihara (2015), Distinguishing time-delayed causal interactions using convergent cross mapping, Sci. Rep. 5 . https://doi.org/10.1038/srep14750. Yu, C., B. Gao, R. Muñoz-Carpena, Y. Tian, L. Wu, and O. Perez-Ovilla (2011), A laboratory study of colloid and solute transport in surface runoff on saturated soil, J. Hydrol. 402 (1), 159-164. https://doi.org/https://doi.org/10.1016/j.jhydrol.2011.03.011. Zhang, D., W. S. Wang, S. Y. Yu, S. Q. Liang, and Q. F. Hu (2021a), Assessment of the Contributions of Climate Change and Human Activities to Runoff Variation: Case Study in Four Subregions of the Jinsha River Basin, China, J. Hydro. Eng. 26 (9). https://doi.org/10.1061/(asce)he.1943-5584.0002119. Zhang, L., E. Moges, J. W. Kirchner, E. Coda, T. C. Liu, A. S. Wymore, Z. X. Xu, and L. G. Larsen (2021b), CHOSEN: A synthesis of hydrometeorological data from intensively monitored catchments and comparative analysis of hydrologic extremes, Hydrol. Process. 35 (11). https://doi.org/10.1002/hyp.14429. Zhao, Y. Y., T. J. Zhu, Z. Q. Zhou, H. J. Cai, and Z. D. Cao (2023a), Detecting nonlinear information about drought propagation time and rate with nonlinear dynamic system and chaos theory, J. Hydrol. 623 . https://doi.org/10.1016/j.jhydrol.2023.129810. Zhao, Y. Y., Y. G. Zou, E. Z. Ma, Z. Q. Zhou, Y. Q. Feng, Z. D. Cao, H. J. Cai, C. Li, and Y. H. Yan (2023b), Can groundwater storage in turn affect the cryospheric variables? A new perspective from nonlinear dynamic causality detection, J. Hydrol. 624 , 14. https://doi.org/10.1016/j.jhydrol.2023.129910. Zounemat-Kermani, M., O. Batelaan, M. Fadaee, and R. Hinkelmann (2021), Ensemble machine learning paradigms in hydrology: A review, J. Hydrol. 598 , 13. https://doi.org/10.1016/j.jhydrol.2021.126266. Additional Declarations No competing interests reported. Cite Share Download PDF Status: Published Journal Publication published 16 Apr, 2025 Read the published version in Stochastic Environmental Research and Risk Assessment → Version 1 posted Editorial decision: Revision requested 23 Sep, 2024 Reviews received at journal 26 Aug, 2024 Reviews received at journal 14 Aug, 2024 Reviewers agreed at journal 30 Jul, 2024 Reviewers agreed at journal 08 Jul, 2024 Reviewers invited by journal 07 Jul, 2024 Editor assigned by journal 28 Jun, 2024 Submission checks completed at journal 26 Jun, 2024 First submitted to journal 26 Jun, 2024 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4643196","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":326343101,"identity":"3a4b94e8-f4a9-4b22-8cb6-4aeca778dc36","order_by":0,"name":"Hanxu Liang","email":"","orcid":"","institution":"Wuhan University","correspondingAuthor":false,"prefix":"","firstName":"Hanxu","middleName":"","lastName":"Liang","suffix":""},{"id":326343103,"identity":"9652e72e-62d3-4789-8d4d-65c643384463","order_by":1,"name":"Wensheng Wang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA00lEQVRIie3LvQrCMBSG4UihXSJd42CvIRAIQgdvpVnaJXZ2cMjUtav3UXAOBHRpcc0qgosOutefiOAY4yaYdzoHvgcAn+93u8P4dQTOJEhG4ltCsHQleNOpI+xD1mzVAYF5ykTUSTtpyzwdVpCtdE4RaAsmYJlZCZWckplAhsAQDSrFBILYTrYnSniPWVO3htxciOZkz8OMYMANEQ5kqk80uFYyQTonk2xdkApyOxnVnFyWvYRxrXb6vEjHddTaiSlE7zN7vp/2puDsMPL5fL5/7gFixEP8cF4wkAAAAABJRU5ErkJggg==","orcid":"","institution":"Sichuan University","correspondingAuthor":true,"prefix":"","firstName":"Wensheng","middleName":"","lastName":"Wang","suffix":""},{"id":326343104,"identity":"1ad59eb0-097e-408b-b08f-fb5602732c3a","order_by":2,"name":"Bin Chen","email":"","orcid":"","institution":"Sichuan Province ABA Hydrology and Water Resources Survey Center","correspondingAuthor":false,"prefix":"","firstName":"Bin","middleName":"","lastName":"Chen","suffix":""},{"id":326343107,"identity":"811825c7-69a7-4824-af85-c177d8c5c3e2","order_by":3,"name":"Li Guo","email":"","orcid":"","institution":"Sichuan University","correspondingAuthor":false,"prefix":"","firstName":"Li","middleName":"","lastName":"Guo","suffix":""},{"id":326343109,"identity":"b3fb59a2-f6e9-4f36-8c92-1e0fe49bff97","order_by":4,"name":"Hu Liu","email":"","orcid":"","institution":"Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Hu","middleName":"","lastName":"Liu","suffix":""},{"id":326343110,"identity":"edd1d7a7-1a2a-4993-9bcb-d346be428394","order_by":5,"name":"Siyi Yu","email":"","orcid":"","institution":"Sichuan University","correspondingAuthor":false,"prefix":"","firstName":"Siyi","middleName":"","lastName":"Yu","suffix":""},{"id":326343113,"identity":"e1f2d4bb-f3ad-49e4-8284-d8e893487a68","order_by":6,"name":"Dan Zhang","email":"","orcid":"","institution":"Changjiang Survey, Design and Research Co., Ltd","correspondingAuthor":false,"prefix":"","firstName":"Dan","middleName":"","lastName":"Zhang","suffix":""}],"badges":[],"createdAt":"2024-06-26 13:52:24","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4643196/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4643196/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1007/s00477-025-02977-3","type":"published","date":"2025-04-16T15:57:25+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":60708202,"identity":"2e4fb24b-82f2-42e6-86de-6c4f49d77ceb","added_by":"auto","created_at":"2024-07-19 19:44:06","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":443020,"visible":true,"origin":"","legend":"\u003cp\u003eOverview of the four causal inference methods. (a) Cross-correlation function. \u0026nbsp;(b) Convergent cross mapping. (c) Transfer entropy. (d) Causal network learning algorithms. The black dots represent the state of variables to be detected; The black dashed box denotes the past of the cause variables \u003cem\u003eX\u003c/em\u003e; The blue solid box denotes the conditional object in conditional independence tests. \u003cem\u003eM\u003c/em\u003e is the original manifold constructed by the trajectories of\u003cem\u003e X\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e\u003cem\u003e,\u003c/em\u003e \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e, and \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e in the phase space; \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ex\u003c/em\u003e\u003c/sub\u003e and \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ey\u003c/em\u003e\u003c/sub\u003e denote the shadow manifold constructed by the current and lagged states of variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e, respectively\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/7295f5034952ec79ace57d1b.png"},{"id":60707752,"identity":"2203ac35-6c66-495a-8979-603b8476d729","added_by":"auto","created_at":"2024-07-19 19:36:06","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":166498,"visible":true,"origin":"","legend":"\u003cp\u003eStructures and causal networks of the conceptual hydrological model and distributions of the generated data. (a) Model Structure 1: rainfall-runoff structure. (b) Model Structure 2: snowmelt-runoff structure.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/2fae4183458355a65bebf65c.png"},{"id":60706274,"identity":"bdffb9ff-1273-440f-92b2-856be4eeac7c","added_by":"auto","created_at":"2024-07-19 19:28:06","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":343668,"visible":true,"origin":"","legend":"\u003cp\u003eThe causal structures of the hydrological model identified by different causal inference methods. A causal link exists only if it is identified by more than half of the 100 simulations.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/9c2c9c20332576db9b47bb75.png"},{"id":60706282,"identity":"759b6777-ffc9-4945-afd6-f40272bc67b9","added_by":"auto","created_at":"2024-07-19 19:28:07","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":289360,"visible":true,"origin":"","legend":"\u003cp\u003eThe impact of sample length on the performance of different causal methods. Dashed lines denote the situations where direct and indirect causality are not distinguished.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/d6a1a35a2d8ff16462d96b35.png"},{"id":60706284,"identity":"3c450d04-2ee2-48b2-a67c-da24cebfecf9","added_by":"auto","created_at":"2024-07-19 19:28:07","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":241594,"visible":true,"origin":"","legend":"\u003cp\u003eThe impact of noise on the performance of different causal methods. Dashed lines denote the situations where direct and indirect causality are not distinguished.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/90ee7acfaf370b0d9ded3361.png"},{"id":60706279,"identity":"a192c552-bde6-41f3-adc9-ef7ec571e621","added_by":"auto","created_at":"2024-07-19 19:28:06","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":587315,"visible":true,"origin":"","legend":"\u003cp\u003eThe impact of missing data on the performance of different causal methods. (a)-(d) represent the synchronous sparsity; (e)-(h) represent the asynchronous sparsity. Dashed lines denote the situations where direct and indirect causality are not distinguished.\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/9c0ed27b6e419d7956e2afa5.png"},{"id":60706283,"identity":"747e391f-1e2d-4fb2-93e8-46af917087ad","added_by":"auto","created_at":"2024-07-19 19:28:07","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":869403,"visible":true,"origin":"","legend":"\u003cp\u003eLocation and data for the Shale Hills Catchment. (a) geographic location of the Shale Hills Catchment. (b) and (c) show the min-max normalized values of the time series during snow-free (May to October) and snow-cover (November to April) periods, respectively.\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/ebf37070b9d50450a90fd37c.png"},{"id":60706280,"identity":"112e0b22-a264-4c2e-9125-63bf704e7bbc","added_by":"auto","created_at":"2024-07-19 19:28:06","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":479949,"visible":true,"origin":"","legend":"\u003cp\u003eThe causal structures for the Shale Hills Catchment during the snow-free period. The curve and straight line represent time-lagged and contemporaneous causal links, respectively.\u003c/p\u003e","description":"","filename":"8.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/e792f0819eba52c272e6b937.png"},{"id":60706277,"identity":"744f7eb9-15c9-4a5c-81d1-80b7861ba455","added_by":"auto","created_at":"2024-07-19 19:28:06","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":287853,"visible":true,"origin":"","legend":"\u003cp\u003eThe causal structures for the Shale Hills Catchment during the snow-cover period. The curve and straight line represent time-lagged and contemporaneous causal links, respectively.\u003c/p\u003e","description":"","filename":"9.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/f6f44ba9e19308f93ff61c15.png"},{"id":60706278,"identity":"deaf4100-7b3f-4ac3-b424-11826516bcb9","added_by":"auto","created_at":"2024-07-19 19:28:06","extension":"png","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":757414,"visible":true,"origin":"","legend":"\u003cp\u003eLocation and data for the Chuosijia River Basin. (a) geographic location of the basin. (b) and (c) show the min-max normalized values of the time series during snow-free (June to October) and snow-cover (November to May) periods, respectively.\u003c/p\u003e","description":"","filename":"10.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/4c6a0a2868dd76eded0ea5dc.png"},{"id":60706285,"identity":"4def9af7-3936-406e-b323-b897aa330ef1","added_by":"auto","created_at":"2024-07-19 19:28:07","extension":"png","order_by":11,"title":"Figure 11","display":"","copyAsset":false,"role":"figure","size":450876,"visible":true,"origin":"","legend":"\u003cp\u003eThe causal structures for the Chuosijia River Basin during the snow-free period. The curve and straight line represent time-lagged and contemporaneous causal links, respectively.\u003c/p\u003e","description":"","filename":"11.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/16854b5f99071d6a7c21b475.png"},{"id":60706281,"identity":"abf1e91a-68c7-495d-a312-0f8d65e4c6e7","added_by":"auto","created_at":"2024-07-19 19:28:07","extension":"png","order_by":12,"title":"Figure 12","display":"","copyAsset":false,"role":"figure","size":612872,"visible":true,"origin":"","legend":"\u003cp\u003eThe causal structures for the Chuosijia River Basin during the snow-cover period. The curve and straight line represent time-lagged and contemporaneous causal links, respectively.\u003c/p\u003e","description":"","filename":"12.png","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/9ab6d21fd13f6731f2f1f821.png"},{"id":81050818,"identity":"25224b52-5379-4895-88fc-82d2170579ac","added_by":"auto","created_at":"2025-04-21 16:05:35","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":8072872,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4643196/v1/c90c897b-c7db-4da5-a61e-87df7c70de0e.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Inferring causal associations in hydrological systems: A comparison of methods","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eThe hydrological systems, which encompass the interaction of multiple variables (precipitation, runoff, evaporation, etc.) across spatiotemporal scales, are complex giant systems with both deterministic and stochastic dynamics (Sang et al., \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Many hydrological questions are inherently causal, which aim to understand whether and how the cause variables impact the effect variables. For instance, how to identify and quantify interactions among variables in the hydrologic cycle (Good et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Kleidon and Renner, \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2013\u003c/span\u003e); how to select appropriate influence factors for robust hydrologic predictions (Apaydin and Sibtain, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Chen et al., \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2011\u003c/span\u003e); how to determine the direct and indirect factors in the analysis of streamflow or hydrologic extremes (i.e., flood and drought) under the changing environment (Liang et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Zhang et al., \u003cspan citationid=\"CR84\" class=\"CitationRef\"\u003e2021a\u003c/span\u003e; Peng et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). It is widely acknowledged that correlation does not imply causation (Altman and Krzywinski, \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). Associations can arise between variables in the presence (i.e., \u003cem\u003eX\u003c/em\u003e causes \u003cem\u003eY\u003c/em\u003e) and absence (i.e., they have a common cause) of a causal relationship, as stated in Reichenbach\u0026rsquo;s common cause principle (Reichenbach, \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e1956\u003c/span\u003e). Relying solely on correlation-based methods may not fully encapsulate the intrinsic causal mechanisms in hydrological systems, especially within the complex spatiotemporal fabric characterized by numerous variables. Besides, Bl\u0026ouml;schl et al. (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2019\u003c/span\u003e) recently stated in the 23 unsolved problems for hydrologic sciences that\u0026ldquo;\u003cem\u003eQuestions remain focused on process-based understanding of hydrological variability and causality at all space and time scales\u003c/em\u003e\u0026rdquo;, which highlights the urgent need for the causal perspective and framework to promote a better understanding of complex hydrological systems.\u003c/p\u003e \u003cp\u003eGenerally, replicated interventional experiments on variables provide an explicit way to understand causal processes. However, it is often unfeasible or unethical in the real world, such as operating and intervening in regional or global hydrological variables. Additionally, with the development of monitoring technology, several experimental catchments have been established in recent years (Zhang et al., \u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e2021b\u003c/span\u003e). However, the construction and maintenance are costly, and due to the spatial heterogeneity and scale problems of hydrological systems (Bergstr\u0026ouml;m and Graham, \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e1998\u003c/span\u003e; Bl\u0026ouml;schl et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), the extensibility and upscaling of discovered mechanisms are also challenging. Moreover, data-informed computer simulations may provide alternative randomized controlled experiments, while these are often time-consuming and inaccurate, and demand massive expert knowledge, which in turn may impose strong mechanistic assumptions on the system (Runge et al., \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2019a\u003c/span\u003e). Fortunately, with the development of observational techniques and data sciences, the amount of available time series data in hydrologic sciences has been increasing over the years, which comes from satellite remote sensing observations, station-based or field site measurements, earth system modeling, and reanalysis products (Li et al., \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2023a\u003c/span\u003e). Such data repositories, along with growing computational efficiency, open up an alternative avenue to use data-driven methods without the intervention of hydrological systems, namely observational causal inference.\u003c/p\u003e \u003cp\u003eOver the past decades, with the progress of statistics and computer science, methods for causal inference from observational time series have been developed rapidly. The most classic work can be traced back to the Granger causality (GC; Granger, \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e1969\u003c/span\u003e), which approves the causal impact from the variable \u003cem\u003eX\u003c/em\u003e to \u003cem\u003eY\u003c/em\u003e if the past of \u003cem\u003eX\u003c/em\u003e contributes to the prediction of the future of \u003cem\u003eY.\u003c/em\u003e However, the validity of the GC framework remains controversial since it only reveals the statistical associations via linear vector autoregression (Shojaie and Fox, \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Later, the nonlinear extension of the original GC framework was developed by integrating information theory, that is, the transfer entropy (TE; Schreiber, \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2000\u003c/span\u003e). Except for the GC framework, some complementary perspectives emerged, from the nonlinear dynamics system based on the reconstruction of attractors, known as the convergent cross mapping (CCM; (Sugihara et al., \u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e2012\u003c/span\u003e)). More recently, the graph-based causal network learning algorithms, such as the Peter\u0026ndash;Clark momentary conditional independence (PCMCI) proposed by Runge et al. (\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e2019b\u003c/span\u003e), become increasingly popular in the causal community. Unlike the data-driven machine learning methods that primarily concentrate on classification or prediction (Zounemat-Kermani et al., \u003cspan citationid=\"CR88\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), the purpose of all causal methods is to discover and quantify the causal interactions of the underlying system based on the measured time series (Runge et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Since explainable machine learning has attracted great interest of researchers recently (Angelov et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2021\u003c/span\u003e), it is still challenging for machine learning, even deep learning models to extract complex mechanisms from this black box (Rudin, \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Thus, causal inference methods are significant for supplementing predictive machine learning to enhance our theoretical cognition of underlying systems (Reichstein et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Over the years, causal inference methods have shown huge potential in numerous fields across Earth system science, for example, the regularity discovery, process understanding, hypothesis validation, and physical model improvement (Su et al., \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; (Faybishenko, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2017\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eUnfortunately, in hydrology, traditional correlation and regression methods are still the most commonly used tools (Massei et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Xu et al., \u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Yu et al., \u003cspan citationid=\"CR83\" class=\"CitationRef\"\u003e2011\u003c/span\u003e), and the application of modern causal methods remains incipient, such as the feedback between rainfall and soil moisture (Luo et al., \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e2018\u003c/span\u003e), the interactions between groundwater and streamflow (Bonotto et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Zhao et al., \u003cspan citationid=\"CR87\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e), and the hydrological connectivity in the karstic areas (Delforge et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Li et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e). The modern causal inference methods do have huge potential in hydrological systems if their underlying assumptions and methodological challenges have been fully considered (Runge et al., \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2019a\u003c/span\u003e). Since limited exploratory studies on the applications of causal inference methods in hydrology have been conducted recently, there still remains much confusion about how and when to apply them, and how to interpret their results, thus requiring a comprehensive comparison and evaluation between them to guide the further popularization in hydrological sciences.\u003c/p\u003e \u003cp\u003eIn this study, four popular methods in the causal inference community, namely cross-correlation function (CCF), convergent cross mapping (CCM), transfer entropy (TE), and a causal network learning algorithm (PCMCI+) were selected, with a detailed explanation of their basic principles and underlying assumptions. Then, the performances of the methods were systematically evaluated in both synthetic and real hydrological systems, thus revealing their potential and limitations for application in hydrology. The remainder of this paper is organized as follows. Section \u003cspan refid=\"Sec2\" class=\"InternalRef\"\u003e2\u003c/span\u003e provides a comprehensive description of the four causal inference methods in our comparison study. Section \u003cspan refid=\"Sec7\" class=\"InternalRef\"\u003e3\u003c/span\u003e successively illustrates the setup of the hydrological model structure, evaluation metrics, and experimental results in the synthetic case, which includes the large sample test and sensitivity analysis. Section \u003cspan refid=\"Sec11\" class=\"InternalRef\"\u003e4\u003c/span\u003e presents the performances of the causal inference methods in two real case studies. Further discussion and conclusions are illustrated in Sections 5 and \u003cspan refid=\"Sec22\" class=\"InternalRef\"\u003e6\u003c/span\u003e, respectively.\u003c/p\u003e"},{"header":"2. Causal inference methods","content":"\u003cp\u003eFour causal inference methods, including the CCF, CCM, TE, and PCMCI+, which are popular and representative in the causal inference community (Runge et al., \u003cspan\u003e2019a\u003c/span\u003e; Runge et al., \u003cspan\u003e2023\u003c/span\u003e), were chosen for our research. The basic characteristics of these methods are listed in Table \u003cspan\u003e1\u003c/span\u003e.\u003c/p\u003e\n\u003cdiv\u003e\n \u003ctable id=\"Tab1\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv\u003eTable 1\u003c/div\u003e\n \u003cdiv\u003e\n \u003cp\u003eSummary of the characteristics of the four causal inference methods adopted in this research.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003ccolgroup cols=\"5\"\u003e\u003c/colgroup\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCross-correlation function\u003c/p\u003e\n \u003cp\u003e(CCF)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eConvergent cross mapping\u003c/p\u003e\n \u003cp\u003e(CCM)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eTransfer\u003c/p\u003e\n \u003cp\u003eentropy\u003c/p\u003e\n \u003cp\u003e(TE)\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCausal network learning algorithm\u003c/p\u003e\n \u003cp\u003e(PCMCI+)\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eLinear/nonlinear\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003elinear\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003enonlinear\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003enonlinear\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003elinear/nonlinear\u003csup\u003ea\u003c/sup\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eType of system\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estochastic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003edeterministic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estochastic\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estochastic\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eDetection of contemporaneous\u003c/p\u003e\n \u003cp\u003ecausal links\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eNo\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eYes/No\u003csup\u003eb\u003c/sup\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eNo\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eYes\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eOther assumptions\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e/\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003enonlinear dynamical system\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e/\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ecausal sufficiency, causal Markov, causal faithfulness\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eReferences\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e(Cryer and Chan, \u003cspan\u003e2008\u003c/span\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e(Sugihara et al., \u003cspan\u003e2012\u003c/span\u003e);\u003c/p\u003e\n \u003cp\u003e(Ye et al., \u003cspan\u003e2015\u003c/span\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e(Schreiber, \u003cspan\u003e2000\u003c/span\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e(Runge, \u003cspan\u003e2020\u003c/span\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003ctfoot\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"5\"\u003eNote: \u003csup\u003ea\u003c/sup\u003eDepends on the algorithms for conditional independence detection. \u003csup\u003eb\u003c/sup\u003eDepends on whether the principle of priority of the cause is adopted.\u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tfoot\u003e\n \u003c/table\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec3\"\u003e\n \u003ch2\u003e2.1 Cross-correlation function\u003c/h2\u003e\n \u003cp\u003eFor a driving variable \u003cem\u003e𝑋\u003c/em\u003e and a response variable \u003cem\u003e𝑌\u003c/em\u003e, causality can be identified through the cross-correlation function (CCF) with the principle of priority of the cause. As shown in Fig. \u003cspan\u003e1\u003c/span\u003ea, the variable \u003cem\u003e𝑋\u003c/em\u003e\u003csub\u003e𝑡\u0026minus;d\u003c/sub\u003e is considered to be the cause of variable \u003cem\u003eY\u003c/em\u003e\u003csub\u003e𝑡\u003c/sub\u003e if the Pearson\u0026rsquo;s correlation coefficient \u0026rho; between \u003cem\u003eX\u003c/em\u003e\u003csub\u003et\u0026thinsp;\u0026minus;\u0026thinsp;d\u003c/sub\u003e and \u003cem\u003eY\u003c/em\u003e\u003csub\u003et\u003c/sub\u003e is significant during their common time domain for at least one value of \u003cem\u003ed\u003c/em\u003e up to \u003cem\u003ed\u003c/em\u003e\u003csub\u003emax\u003c/sub\u003e:\u003c/p\u003e\n \u003cdiv id=\"Equ1\"\u003e\n \u003cdiv id=\"FileID_Equ1\" name=\"EquationSource\"\u003e\u003cimg src=\"https://myfiles.space/user_files/122228_c8a1650c59388082/122228_custom_files/img1721320823.png\"\u003e\u003cbr\u003e\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cem\u003eCov\u003c/em\u003e(\u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u0026minus;d\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e) represents the cross-covariance at lag time \u003cem\u003ed\u003c/em\u003e; \u003cem\u003e\u0026sigma;\u003c/em\u003e\u003csub\u003e\u003cem\u003ex\u003c/em\u003e\u003c/sub\u003e and \u003cem\u003e\u0026sigma;\u003c/em\u003e\u003csub\u003e\u003cem\u003ey\u003c/em\u003e\u003c/sub\u003e denote the standard deviations of time series \u003cem\u003ex\u003c/em\u003e and \u003cem\u003ey\u003c/em\u003e, respectively. The significance of the hypothesis (i.e., 𝜌 is different from zero) is usually evaluated using the Student\u0026rsquo;s t-test (Student, \u003cspan\u003e1908\u003c/span\u003e) with a \u003cem\u003ep\u003c/em\u003e-value. Given a significance level 𝛼, significant relationships are considered if the \u003cem\u003ep\u003c/em\u003e-value is lower than 𝛼.\u003c/p\u003e\n \u003cp\u003eIn this study, significant correlations identified by the CCF method are considered as a way to reveal potentially causal links between two time series at a causal delay. However, according to Reichenbach\u0026rsquo;s common cause principle (Reichenbach, \u003cspan\u003e1956\u003c/span\u003e), this linear correlation may not imply causation, owing to the interference of confounders, autocorrelation of variables, and nonlinear dynamical associations (Dean and Dunsmuir, \u003cspan\u003e2016\u003c/span\u003e). Despite these inherent limitations as a bivariate linear method, the CCF method has still been widely applied in many disciplines including hydrology due to its simplicity, efficiency, and linear interpretability (Massei et al., \u003cspan\u003e2006\u003c/span\u003e; Xu et al., \u003cspan\u003e2022\u003c/span\u003e; Yu et al., \u003cspan\u003e2011\u003c/span\u003e).\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec4\"\u003e\n \u003ch2\u003e2.2 Convergent cross mapping\u003c/h2\u003e\n \u003cp\u003eConvergent cross mapping (CCM), proposed by Sugihara et al. (\u003cspan\u003e2012\u003c/span\u003e), aims to infer causality between two time series in nonlinear dynamical systems using the theory of time delay embedding (Takens, \u003cspan\u003e1981\u003c/span\u003e). It detects topological similarities between reconstructed attractor manifolds (i.e., consistent behaviors when revisiting similar states). The basic assumption of the CCM is that variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e belong to the same nonlinear dynamical system, and if \u003cem\u003eX\u003c/em\u003e causes \u003cem\u003eY\u003c/em\u003e, the cause variable \u003cem\u003eX\u003c/em\u003e leaves informative features in the affected variable \u003cem\u003eY\u003c/em\u003e, and thus \u003cem\u003eX\u003c/em\u003e can be estimated from its features in \u003cem\u003eY\u003c/em\u003e (Fig. \u003cspan\u003e1\u003c/span\u003eb).\u003c/p\u003e\n \u003cp\u003eSpecifically, to identify the causation between variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e, the first stage is to reconstruct the shadow manifold by the lagged time series of each variable. For instance, the trajectories \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ex\u003c/em\u003e\u003c/sub\u003e={\u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u0026minus;\u0026tau;\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;2\u003cem\u003e\u0026tau;\u003c/em\u003e\u003c/sub\u003e} and \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ey\u003c/em\u003e\u003c/sub\u003e={\u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u0026minus;\u0026tau;\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u0026minus;\u003c/em\u003e2\u003cem\u003e\u0026tau;\u003c/em\u003e\u003c/sub\u003e} are the shadow manifold of variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e, respectively, and \u003cem\u003e\u0026tau;\u003c/em\u003e denotes the time delay. Secondly, for a single point in \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ey\u003c/em\u003e\u003c/sub\u003e, \u003cem\u003ek\u003c/em\u003e (=\u0026thinsp;\u003cem\u003eE\u003c/em\u003e\u0026thinsp;+\u0026thinsp;1) nearest neighbors \u003cspan\u003e\u003cspan\u003e\\({Y_{{t_m}}}\\)\u003c/span\u003e\u003c/span\u003e(\u003cem\u003em\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1, 2,\u0026hellip;\u003cem\u003ek\u003c/em\u003e) are identified based on their Euclidean distance to \u003cem\u003eY\u003c/em\u003e, where \u003cem\u003eE\u003c/em\u003e is the embedding dimension. Next, time indices of these identified nearest neighbors are utilized to map corresponding points \u003cspan\u003e\u003cspan\u003e\\({X_{{t_m}}}\\)\u003c/span\u003e\u003c/span\u003ein \u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003ex\u003c/em\u003e\u003c/sub\u003e, and their weighted average values considering the Euclidean distances are calculated to estimate \u003cem\u003eX\u003c/em\u003e, shown as follows:\u003c/p\u003e\n \u003cdiv id=\"Equ2\"\u003e\n \u003cdiv id=\"FileID_Equ2\" name=\"EquationSource\"\u003e$${\\mathop X\\limits^{ \\wedge } _t}|{M_y}=\\sum\\limits_{{m=1}}^{k} {{\\omega _m}} {X_{{t_m}}}$$\u003c/div\u003e\n \u003cdiv\u003e2\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ3\"\u003e\n \u003cdiv id=\"FileID_Equ3\" name=\"EquationSource\"\u003e$${\\omega _m}=\\frac{{{u_m}}}{{\\sum\\limits_{{n=1}}^{k} {{u_n}} }}$$\u003c/div\u003e\n \u003cdiv\u003e3\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ4\"\u003e\n \u003cdiv id=\"FileID_Equ4\" name=\"EquationSource\"\u003e$${u_m}={e^{ - \\frac{{d({Y_t},{Y_{{t_m}}})}}{{d({Y_t},{Y_{{t_1}}})}}}}$$\u003c/div\u003e\n \u003cdiv\u003e4\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cspan\u003e\u003cspan\u003e\\({\\mathop X\\limits^{ \\wedge } _t}\\)\u003c/span\u003e\u003c/span\u003erepresents the estimated values of \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e from \u003cem\u003eM\u003c/em\u003e\u003csub\u003ey\u003c/sub\u003e, and \u003cem\u003ed\u003c/em\u003e(\u003cem\u003eX\u003c/em\u003e, \u003cem\u003eY\u003c/em\u003e) denotes the Euclidean distance between vectors \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e. Finally, the cross-mapping skill, represented by Pearson\u0026apos;s correlation coefficient 𝜌 between observed \u003cem\u003eX\u003c/em\u003e\u003csub\u003et\u003c/sub\u003e and the estimated \u003cspan\u003e\u003cspan\u003e\\({\\mathop X\\limits^{ \\wedge } _t}\\)\u003c/span\u003e\u003c/span\u003e, quantifies the strength of causality from variable \u003cem\u003eX\u003c/em\u003e to \u003cem\u003eY\u003c/em\u003e. Besides, 𝜌 will increase from 0 to a convergent value as the cross-mapping process undergo all points of the shadow manifold seriatim. Referred to Ombadi et al. (\u003cspan\u003e2020\u003c/span\u003e), the significance threshold of 𝜌 was set as 0.3.\u003c/p\u003e\n \u003cp\u003eGiven that the CCM method is based on the chaos theory, it should be applied under the strict assumptions of the deterministic mathematical models, namely, low-dimensional systems with infinite length, no noise, and acyclic and non-intermittent relations (Khatibi et al., \u003cspan\u003e2012\u003c/span\u003e; Sivakumar, \u003cspan\u003e2004\u003c/span\u003e). Additionally, according to Sugihara et al. (\u003cspan\u003e2017\u003c/span\u003e), convergence can be achieved when variables are subject to synchrony, a well-known phenomenon resulting from strong coupling or dynamic resonance, which might lead to spurious bidirectional causal relationships, especially for hydrological variables (Ombadi et al., \u003cspan\u003e2020\u003c/span\u003e; Zhao et al., \u003cspan\u003e2023a\u003c/span\u003e). Therefore, in this study, the time-lagged CCM (Ye et al., \u003cspan\u003e2015\u003c/span\u003e), which aims to overcome the synchrony with the asymmetric patterns of time dependencies (i.e., the principle of priority of the cause), was adopted for our research. Distinguished from the traditional CCM method that uses \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e to cross map \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e, the time-lagged CCM method adopts \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et+d\u003c/em\u003e\u003c/sub\u003e instead. Note that \u003cem\u003ed\u003c/em\u003e is the lag time of cross-mapping, ranging from 0 to \u003cem\u003ed\u003c/em\u003e\u003csub\u003emax\u003c/sub\u003e, which is different from the embedding delay \u003cem\u003e\u0026tau;\u003c/em\u003e in reconstruction of the shadow manifold.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec5\"\u003e\n \u003ch2\u003e2.3 Transfer entropy\u003c/h2\u003e\n \u003cp\u003eTransfer Entropy (TE) is an information-theoretic approach that quantifies the temporally asymmetric transfer of information between two time series (Schreiber, \u003cspan\u003e2000\u003c/span\u003e). It can be deemed as the nonparametric extension of the Granger Causality (Granger, \u003cspan\u003e1969\u003c/span\u003e) and be adopted to infer the causal association from \u003cem\u003eX\u003c/em\u003e to \u003cem\u003eY\u003c/em\u003e if knowing the information about the past of variable \u003cem\u003eX\u003c/em\u003e can reduce the uncertainty of variable \u003cem\u003eY\u003c/em\u003e. Due to its no assumptions on the underlying structure of data and applicability to nonlinear separable systems, the TE method has been widely applied in hydrology (Bennett et al., \u003cspan\u003e2019\u003c/span\u003e; Konapala et al., \u003cspan\u003e2020\u003c/span\u003e; Moges et al., \u003cspan\u003e2022\u003c/span\u003e). To calculate TE, the concept of mutual information is introduced (Kinney and Atwal, \u003cspan\u003e2014\u003c/span\u003e), which can be considered as the overlapping information between two random variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e, shown as follows:\u003c/p\u003e\n \u003cdiv id=\"Equ5\"\u003e\n \u003cdiv id=\"FileID_Equ5\" name=\"EquationSource\"\u003e$$I(X;Y)=H(X)+H(Y) - H(X,Y)$$\u003c/div\u003e\n \u003cdiv\u003e5\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cem\u003eI(X; Y)\u003c/em\u003e denotes the mutual information between variables \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e; \u003cem\u003eH\u003c/em\u003e is the Shannon entropy (Shannon, \u003cspan\u003e1948\u003c/span\u003e), a measure of a variable\u0026rsquo;s average uncertainty; \u003cem\u003eH(X)\u003c/em\u003e and \u003cem\u003eH(X, Y)\u003c/em\u003e represent the univariate and bivariate information entropy, respectively, calculated as:\u003c/p\u003e\n \u003cdiv id=\"Equ6\"\u003e\n \u003cdiv id=\"FileID_Equ6\" name=\"EquationSource\"\u003e$$H(X)= - \\int\\limits_{X} {p(x){{\\log }_2}(p(x))dx}$$\u003c/div\u003e\u003cdiv\u003e6\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equ7\"\u003e\u003cdiv id=\"FileID_Equ7\" name=\"EquationSource\"\u003e$$H(X,Y)= - \\int\\limits_{Y} {\\int\\limits_{X} {p(x,y){{\\log }_2}(p(x,y))dxdy} }$$\u003c/div\u003e\u003cdiv\u003e7\u003c/div\u003e\u003c/div\u003e\u003cp\u003ewhere \u003cem\u003ep\u003c/em\u003e represents the probability distribution of the random variable. Then, TE can be understood as the shared mutual information between the past of the cause variable \u003cem\u003eX\u003c/em\u003e and the present of the affected variable \u003cem\u003eY\u003c/em\u003e, conditioned on the past of \u003cem\u003eY\u003c/em\u003e (Fig. \u003cspan\u003e1\u003c/span\u003ec), calculated as:\u003c/p\u003e\u003cdiv id=\"Equ8\"\u003e\u003cdiv id=\"FileID_Equ8\" name=\"EquationSource\"\u003e$$T{E_{X \\to Y}}=I({Y_t};{X_{t - d}}|{Y_{t - 1}})$$\u003c/div\u003e\n \u003cdiv\u003e8\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ9\"\u003e\n \u003cdiv id=\"FileID_Equ9\" name=\"EquationSource\"\u003e$$I(X;Y|Z)=H(XZ)+H(Y|Z) - H(X,YZ)$$\u003c/div\u003e\n \u003cdiv\u003e9\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ10\"\u003e\n \u003cdiv id=\"FileID_Equ10\" name=\"EquationSource\"\u003e$$H(X|Y)=H(X) - I(X;Y)$$\u003c/div\u003e\n \u003cdiv\u003e10\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cem\u003ed\u003c/em\u003e is the lag time in the transfer of information between \u003cem\u003eX\u003c/em\u003e and \u003cem\u003eY\u003c/em\u003e; \u003cem\u003eI(X;Y|Z)\u003c/em\u003e and \u003cem\u003eH(X|Y)\u003c/em\u003e denote the conditional mutual information and conditional entropy, respectively. Besides, \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u0026minus;\u003c/em\u003e1\u003c/sub\u003e was selected to represent the past of variable \u003cem\u003eY\u003c/em\u003e based on the assumption of the Markov process that the immediate history always provides the most information (Budakoti et al., \u003cspan\u003e2021\u003c/span\u003e; Ruddell and Kumar, \u003cspan\u003e2009\u003c/span\u003e). To estimate the probability distribution of the variables, following Ruddell and Kumar (\u003cspan\u003e2009\u003c/span\u003e), a histogram-based approach was adopted with 11 bins. Considering that the hydrological time series includes several zero values, which may cause deviation in the calculation of entropy, we used the approach proposed by Gong et al. (\u003cspan\u003e2014\u003c/span\u003e) to handle the zero and nonzero values separately in the binning. Additionally, the significance of TE was tested by the Monte Carlo methods and Student\u0026rsquo;s t-test. The TE values were in comparison with a series of 500 randomly shuffled surrogates in which any temporal correlations between the two time series were broken (Moges et al., \u003cspan\u003e2022\u003c/span\u003e). The TE value is considered significant if it exceeds the threshold, defined as the value corresponding to the 95th percentile of the random samples.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec6\"\u003e\n \u003ch2\u003e2.4 Causal network learning algorithms\u003c/h2\u003e\n \u003cp\u003eThe causal network learning algorithms, based on a series of graphical rules that dominate the identification of system causal associations, have been developed for reconstructing large-scale causal graphs with high dimensionality of variables (Runge et al., \u003cspan\u003e2019b\u003c/span\u003e). To infer the structure of the causal graph, three underlying assumptions, namely causal sufficiency, causal Markov, and causal faithfulness are usually required (Runge, \u003cspan\u003e2018\u003c/span\u003e). The causal sufficiency assumes that there are no other unobserved (or latent) variables that directly or indirectly impact any pair of our variable sets (i.e. all variables are directly observed). The other two assumptions allow to relate \u003cem\u003ed-separations\u003c/em\u003e to conditional independence in the graph under this distribution: \u003cem\u003eX d-sep Y\u003c/em\u003e|\u003cem\u003eZ\u003c/em\u003e ⟺ \u003cem\u003eX\u003c/em\u003e\u0026perp;\u003cem\u003eY\u003c/em\u003e|\u003cem\u003eZ\u003c/em\u003e, where causal Markov and causal faithfulness assumptions correspond to the forward and backward arrows, respectively. Two nodes can be deemed as \u003cem\u003ed-separated\u003c/em\u003e given \u003cem\u003eZ\u003c/em\u003e if and only if all paths are blocked given \u003cem\u003eZ\u003c/em\u003e (Pearl, \u003cspan\u003e2009\u003c/span\u003e). The causal faithfulness assumption precludes the case where \u003cem\u003eX\u003c/em\u003e affects \u003cem\u003eY\u003c/em\u003e in two directions with positive and negative effects canceling out.\u003c/p\u003e\n \u003cp\u003eIn this study, the latest causal network learning algorithm, namely Peter-Clark momentary conditional independence plus (PCMCI+), which solves the problems for detecting contemporaneous and time-lagged causal links in autocorrelated high dimensional time series (Runge, \u003cspan\u003e2020\u003c/span\u003e), was adopted. The algorithm can be divided into two stages, namely the skeleton discovery phase and the orientation phase (Fig. \u003cspan\u003e1\u003c/span\u003ed). In the skeleton discovery phase stage, the Peter\u0026ndash;Clark (PC) algorithm begins with a fully connected graph and then tests for the elimination of links between variables iteratively by conditioning sets of increasing cardinality (Spirtes and Glymour, 1991). Specifically, it starts with a graph where all unconditionally (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0) dependent variable pairs are connected with the assumption of the stationarity for causal links, and iteratively tests conditional independence with an increasing number of conditions \u003cem\u003ep\u003c/em\u003e. Lagged links are oriented forward in time based on the principle of the priority of cause, whereas contemporaneous links are left undirected (circle marks at the ends) in this skeleton discovery phase. For instance, \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;1\u003c/sub\u003e and \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e (black nodes) have already been correctly identified as independent in the second iteration step (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1) where the dependence through \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;1\u003c/sub\u003e (blue box) is conditioned out, while we need to condition on two variables to test whether \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;2\u003c/sub\u003e and \u003cem\u003eW\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e are independent (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2). To further eliminate spurious links for all ordered pairs and increase the detection power, the contemporaneous conditions are iterated over by momentary conditional independence (MCI) tests with partial correlation:\u003c/p\u003e\n \u003cp\u003e\u003cspan\u003e\u0026nbsp;\u003c/span\u003e \u003cspan\u003e\u003cspan\u003e\\(MCI:X_{{t - d}}^{i} \\bot X_{t}^{j}|\\mathop P\\limits^{ \\wedge } (X_{t}^{j})\\backslash \\{ X_{{t - d}}^{i}\\} ,\\mathop P\\limits^{ \\wedge } (X_{{t - d}}^{i})\\)\u003c/span\u003e\u003c/span\u003e (11)\u003c/p\u003e\n \u003cp\u003ewhere \u003cspan\u003e\u003cspan\u003e\\(\\widehat{P}\\left({X}_{t}^{j}\\right)\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan\u003e\u003cspan\u003e\\(\\widehat{P}\\left({X}_{t-d}^{i}\\right)\\)\u003c/span\u003e\u003c/span\u003e represent the lagged parents of \u003cspan\u003e\u003cspan\u003e\\({X}_{t}^{j}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan\u003e\u003cspan\u003e\\({X}_{t-d}^{i}\\)\u003c/span\u003e\u003c/span\u003e identified in the previous PC step, respectively; \u0026lsquo;\\\u0026rsquo; means precluding \u003cspan\u003e\u003cspan\u003e\\({X}_{t-d}^{i}\\)\u003c/span\u003e\u003c/span\u003e from \u003cspan\u003e\u003cspan\u003e\\(\\widehat{P}\\left({X}_{t}^{j}\\right)\\)\u003c/span\u003e\u003c/span\u003e. After the iteration through subsets \u003cem\u003eS\u003c/em\u003e \u003cspan\u003e\u003cspan\u003e\\(\\subset\\)\u003c/span\u003e\u003c/span\u003e \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e of contemporaneous links, the spurious adjacencies can be fully removed. The MCI partial correlation values were used to represent the link strength of causality, with the significance level \u003cem\u003e\u0026alpha;\u003c/em\u003e\u003csub\u003epc\u003c/sub\u003e set to 0.05 in the Student\u0026rsquo;s t-test.\u003c/p\u003e\n \u003cp\u003eIn the second orientation phase, the left contemporaneous links can finally be directed by a series of rules. For example, when finding that \u003cem\u003eW\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;1\u003c/sub\u003e and \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e are independent conditional on \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u0026minus;1\u003c/sub\u003e, while not independent conditional on \u003cem\u003eW\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e, the causal links from \u003cem\u003eZ\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e to \u003cem\u003eW\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e can be identified since the other causal directions are not in accordance with the observed conditional independencies. However, there may not exist such rules to distinguish the direction between the \u003cem\u003eX\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e and \u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003et\u003c/em\u003e\u003c/sub\u003e due to the Markov equivalence class (Pearl, \u003cspan\u003e2009\u003c/span\u003e). More details about the PCMCI\u0026thinsp;+\u0026thinsp;algorithm, along with its pseudo code, can be found in Runge (\u003cspan\u003e2020\u003c/span\u003e). To the best of our knowledge, the PCMCI\u0026thinsp;+\u0026thinsp;algorithm has not been applied in hydrological systems hitherto, and this study provides the first evaluation of its potential for further applications.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"3. Application to synthetic case study","content":"\u003cdiv id=\"Sec8\"\u003e\n \u003ch2\u003e3.1 Construction of conceptual hydrological model\u003c/h2\u003e\n \u003cp\u003eTo evaluate the ability of the above four methods for identifying causal associations in hydrological systems, firstly, a conceptual hydrological model was employed to generate the synthetic time series, where the underlying causal associations have already been acquired and can be deemed as true relationships. The EXP-HYDRO model (Patil and Stieglitz, \u003cspan\u003e2014\u003c/span\u003e), which characterizes the complexity of hydrological processes while maintaining as few variables and parameters as possible to ensure the interpretability of the constructed causal network, was chosen as our basic model framework. We simplified and split the original EXP-HYDRO model into two structures, namely the rainfall-runoff structure (Model S1, Fig,2a) and the snowmelt-runoff structure (Model S2, Fig. \u003cspan\u003e2\u003c/span\u003eb).\u003c/p\u003e\n \u003cp\u003eModel S1 includes four variables: effective precipitation \u003cem\u003eP\u003c/em\u003e, soil water storage \u003cem\u003eS\u003c/em\u003e, interflow \u003cem\u003eI\u003c/em\u003e, and runoff \u003cem\u003eQ\u003c/em\u003e. Besides, the model contains three parameters: the maximum soil water storage \u003cem\u003eS\u003c/em\u003e\u003csub\u003emax\u003c/sub\u003e, and nonlinear storage-discharge parameters \u003cem\u003eK\u003c/em\u003e\u003csub\u003e\u003cem\u003es\u003c/em\u003e,\u003c/sub\u003e and \u003cem\u003e\u0026gamma;\u003c/em\u003e. The model structure and corresponding causal network are shown at the top of Fig. \u003cspan\u003e2\u003c/span\u003ea. The only forcing variable \u003cem\u003eP\u003c/em\u003e was generated by a stochastic weather generator (WeaGETS; Chen et al., \u003cspan\u003e2010\u003c/span\u003e), including a Markov chain for occurrence and a gamma distribution for quantity. Specifically, the occurrence of precipitation was generated with a binary variable \u003cspan\u003e\u003cspan\u003e\\(\\widehat {P}\\)\u003c/span\u003e\u003c/span\u003e (i.e., wet or dry days) using the first-order Markov chain (Richardson, \u003cspan\u003e1981\u003c/span\u003e):\u003c/p\u003e\n \u003cp\u003e\u003cspan\u003e\u0026nbsp;\u003c/span\u003e \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{i,j}}(t)=probability({\\widehat {P}_t}=j|{\\widehat {P}_{t - 1}}=i){\\text{ ; }}i,j=0,{\\text{ }}1{\\text{ ; }}t\u0026gt;1\\)\u003c/span\u003e\u003c/span\u003e (12)\u003c/p\u003e\n \u003cp\u003ewhere \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{i,j}}\\)\u003c/span\u003e\u003c/span\u003eis the transition probability from the state \u003cem\u003ei\u003c/em\u003e to \u003cem\u003ej\u003c/em\u003e; the states 1 and 0 represent the wet and dry days, respectively. Considering that precipitation either occurs or not on a given day, \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{0,1}}+{\\widetilde {p}_{0,0}}={\\widetilde {p}_{1,0}}+{\\widetilde {p}_{1,1}}=1\\)\u003c/span\u003e\u003c/span\u003e. In this study, we set \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{0,0}}={\\widetilde {p}_{1,1}}=0.8\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{0,1}}={\\widetilde {p}_{1,0}}=0.2\\)\u003c/span\u003e\u003c/span\u003e. The values of \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{0,0}}\\)\u003c/span\u003e\u003c/span\u003eand \u003cspan\u003e\u003cspan\u003e\\({\\widetilde {p}_{1,1}}\\)\u003c/span\u003e\u003c/span\u003e were set relatively high in order to raise the persistence in the model, which could make some obstacles to the detection of causal associations. Next, the quantity of precipitation was generated by a Gamma distribution:\u003cspan\u003e\u003cspan\u003e\\(\\widehat {P}\\sim Gamma(\\alpha ,\\beta)\\);\u003c/span\u003e\u003c/span\u003e Generally, for precipitation data, the parameter \u003cem\u003e\u0026alpha;\u003c/em\u003e is smaller than 1 while \u003cem\u003e\u0026beta;\u003c/em\u003e has a wide range of possible values (Geng et al., \u003cspan\u003e1986\u003c/span\u003e). In this study, \u003cem\u003e\u0026alpha;\u003c/em\u003e and \u003cem\u003e\u0026beta;\u003c/em\u003e were set as 0.6, and 6, respectively. Besides, adding some noise to the forcing variable helps to approach the real situation and test the anti-interference ability for causal inference methods. Here, \u003cspan\u003e\u003cspan\u003e\\({P_t}={\\widehat {P}_t}+{\\eta _p}\\)\u003c/span\u003e\u003c/span\u003e, and \u003cspan\u003e\u003cspan\u003e\\({\\eta _p}\\)\u003c/span\u003e\u003c/span\u003eis the white Gaussian noise (zero mean and \u003cspan\u003e\u003cspan\u003e\\(\\sigma _{P}^{2}\\)\u003c/span\u003e\u003c/span\u003evariance). The \u003cspan\u003e\u003cspan\u003e\\(\\sigma _{P}^{2}\\)\u003c/span\u003e\u003c/span\u003e can be calculated as the expectation value of \u003cem\u003eP\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e divided by the signal-to-noise ratio (\u003cem\u003eSNR\u003c/em\u003e), as follows: \u003cspan\u003e\u003cspan\u003e\\({\\sigma _P}^{2}=E{\\text{ }}[{P^2}]/SNR\\)\u003c/span\u003e\u003c/span\u003e and the \u003cem\u003eSNR\u003c/em\u003e can be transformed into the unit of dB:\u003cspan\u003e\u003cspan\u003e\\(SNR(dB)=10{\\log _{10}}SNR\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e\n \u003cp\u003eThe other variables \u003cem\u003eS\u003c/em\u003e, \u003cem\u003eI\u003c/em\u003e, and \u003cem\u003eQ\u003c/em\u003e were determined by the equations of water balance and nonlinear storage-discharge relationship, as:\u003c/p\u003e\n \u003cdiv id=\"Equ11\"\u003e\n \u003cdiv id=\"FileID_Equ11\" name=\"EquationSource\"\u003e$$\\frac{{d{S_t}}}{{dt}}={P_{t - 1}} - {I_{t - 1}}$$\u003c/div\u003e\n \u003cdiv\u003e13\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ12\"\u003e\n \u003cdiv id=\"FileID_Equ12\" name=\"EquationSource\"\u003e$${I_t}={K_s} \\times S_{{t - 1}}^{\\gamma }$$\u003c/div\u003e\n \u003cdiv\u003e14\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ13\"\u003e\n \u003cdiv id=\"FileID_Equ13\" name=\"EquationSource\"\u003e$${Q_t}=\\left\\{ \\begin{gathered} {S_{t - 1}}+{P_{t - 1}} - {S_{\\hbox{max} }}{\\text{ }};{\\text{ }}{S_t} \\geqslant {S_{\\hbox{max} }} \\hfill \\\\ {\\text{ }}0{\\text{ ; }}{S_t}\u0026lt;{S_{\\hbox{max} }} \\hfill \\\\ \\end{gathered} \\right.$$\u003c/div\u003e\n \u003cdiv\u003e15\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere \u003cem\u003eK\u003c/em\u003e\u003csub\u003e\u003cem\u003es\u003c/em\u003e\u003c/sub\u003e represents the speed of reduction in water storage, and \u003cem\u003e\u0026gamma;\u003c/em\u003e is an exponent for \u003cem\u003eS.\u003c/em\u003e\u003c/p\u003e\n \u003cp\u003eModel S2 includes four variables: temperature \u003cem\u003eT\u003c/em\u003e, snowmelt \u003cem\u003eM\u003c/em\u003e, soil water storage \u003cem\u003eS\u003c/em\u003e, and interflow \u003cem\u003eI\u003c/em\u003e, and three parameters: the nonlinear storage-discharge parameters \u003cem\u003eK\u003c/em\u003e\u003csub\u003e\u003cem\u003es\u003c/em\u003e,\u003c/sub\u003e and \u003cem\u003e\u0026gamma;\u003c/em\u003e, and snow-melting parameter \u003cem\u003eK\u003c/em\u003e\u003csub\u003e\u003cem\u003emelt\u003c/em\u003e\u003c/sub\u003e. The model structure and corresponding causal network are shown at the top of Fig. \u003cspan\u003e2\u003c/span\u003eb. The only forcing variable \u003cem\u003eT\u003c/em\u003e was generated by a Normal distribution with a mean of zero. We also added some noise to this forcing variable: \u003cspan\u003e\u003cspan\u003e\\({T_t}={\\widehat {T}_t}+{\\eta _T}\\)\u003c/span\u003e\u003c/span\u003e, where \u003cspan\u003e\u003cspan\u003e\\({\\eta _T}\\)\u003c/span\u003e\u003c/span\u003eis the white Gaussian noise (zero mean and \u003cspan\u003e\u003cspan\u003e\\(\\sigma _{T}^{2}\\)\u003c/span\u003e\u003c/span\u003evariance). The calculation of \u003cspan\u003e\u003cspan\u003e\\(\\sigma _{T}^{2}\\)\u003c/span\u003e\u003c/span\u003e is similar to \u003cspan\u003e\u003cspan\u003e\\(\\sigma _{P}^{2}\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e\n \u003cp\u003eThe other variables \u003cem\u003eS\u003c/em\u003e, \u003cem\u003eI\u003c/em\u003e, and \u003cem\u003eM\u003c/em\u003e were determined by the equations of water balance and nonlinear storage-discharge relationship, as:\u003c/p\u003e\n \u003cdiv id=\"Equ14\"\u003e\n \u003cdiv id=\"FileID_Equ14\" name=\"EquationSource\"\u003e$${I_t}={K_s} \\times S_{{t - 1}}^{\\gamma }$$\u003c/div\u003e\n \u003cdiv\u003e16\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ15\"\u003e\n \u003cdiv id=\"FileID_Equ15\" name=\"EquationSource\"\u003e$$\\frac{{d{S_t}}}{{dt}}={M_{t - 1}} - {I_{t - 1}}$$\u003c/div\u003e\n \u003cdiv\u003e17\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ16\"\u003e\n \u003cdiv id=\"FileID_Equ16\" name=\"EquationSource\"\u003e$${M_t}=\\left\\{ \\begin{gathered} {K_{melt}} \\times {T_{t - 1}}{\\text{ }};{\\text{ }}{T_{t - 1}} \\geqslant 0 \\hfill \\\\ {\\text{ }}0{\\text{ ; }}{T_{t - 1}}\u0026lt;0 \\hfill \\\\ \\end{gathered} \\right.$$\u003c/div\u003e\n \u003cdiv\u003e18\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003eThe parameters of the conceptual hydrological model are listed in Table \u003cspan\u003e2\u003c/span\u003e and the corresponding frequency distributions of the generated hydrological variables are shown at the bottom of Fig. \u003cspan\u003e2\u003c/span\u003e. The frequency distributions are estimated with a sample size of 2000 and \u003cem\u003eSNR\u003c/em\u003e(dB) of 40. Since most hydrological series follow a skewed distribution and have considerable zeros values (Gong et al., \u003cspan\u003e2014\u003c/span\u003e), especially for the forcing variable rainfall/snow, the synthetic series match these properties and can represent the real hydrological system to some extent.\u003c/p\u003e\n \u003cp\u003eBesides, for each model structure, 100 datasets were generated to ensure the robustness of the assessment. In the sensitivity evaluation of causal methods, sample size varies from 100 to 2000; \u003cem\u003eSNR\u003c/em\u003e(dB) from 2 to 40; data missing rate from 10\u0026ndash;60% with two scenarios, namely synchronous and asynchronous sparsity. In the former scenario, all variables are missing at the same time indices that are generated randomly, while in the latter, the time indices vary from the variables. The missing values are complemented by the arithmetic mean of observed values of the same variable (Gao et al., \u003cspan\u003e2018\u003c/span\u003e).\u003c/p\u003e\n \u003cdiv\u003e \u0026nbsp;\u003ctable id=\"Tab2\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv\u003eTable 2\u003c/div\u003e\n \u003cdiv\u003e\n \u003cp\u003eThe parameters of the conceptual hydrological model.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eParameter\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eValue\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003emaximum soil water storage (\u003cem\u003eS\u003c/em\u003e\u003csub\u003emax\u003c/sub\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e20\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eModel S1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estorage-discharge parameter 1 (\u003cem\u003eK\u003c/em\u003e\u003csub\u003es\u003c/sub\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.01\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estorage-discharge parameter 2 (\u003cem\u003e\u0026gamma;\u003c/em\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e1.5\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003esnow-melting parameter (\u003cem\u003eK\u003c/em\u003e\u003csub\u003emelt\u003c/sub\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e2\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eModel S2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estorage-discharge parameter 1 (\u003cem\u003eK\u003c/em\u003e\u003csub\u003es\u003c/sub\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.5\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003estorage-discharge parameter 2 (\u003cem\u003e\u0026gamma;\u003c/em\u003e)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.6\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n \u003cp\u003eIn this study, three metrics, namely True Positives Rate (TPR), False Positives Rate (FPR), and Accuracy Rate (AR), were adopted to evaluate the performance of causal inference methods:\u003c/p\u003e\n \u003cdiv id=\"Equ17\"\u003e\n \u003cdiv id=\"FileID_Equ17\" name=\"EquationSource\"\u003e$$TPR=\\frac{{TP}}{{TP+FN}}$$\u003c/div\u003e\n \u003cdiv\u003e19\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ18\"\u003e\n \u003cdiv id=\"FileID_Equ18\" name=\"EquationSource\"\u003e$$FPR=\\frac{{FP}}{{FP+TN}}$$\u003c/div\u003e\n \u003cdiv\u003e20\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Equ19\"\u003e\n \u003cdiv id=\"FileID_Equ19\" name=\"EquationSource\"\u003e$$AR=\\frac{{TP+TN}}{{TP+FP+TN+FN}}$$\u003c/div\u003e\n \u003cdiv\u003e21\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003ewhere TP and FN denote the true positive and false negative, respectively, which means that the causal link generated by the model is correctly and incorrectly identified, respectively; TN and FP are true negative and false positive, respectively, which means that the non-causal link is correctly and incorrectly identified, respectively. The lower values of FPR and the higher values of TPR and AR indicate the better performance of the tested causal method.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec9\"\u003e\n \u003ch2\u003e3.2 Performance in large sample tests\u003c/h2\u003e\n \u003cp\u003eIn this section, 100 datasets were generated with a sample length of 2000 to evaluate the performance of different causal methods. Figure \u003cspan\u003e3\u003c/span\u003e presents the causal structures of the hydrological model identified by four causal methods. The results represent the average behavior of the experiment, i.e., a causal link exists only if it is identified by more than half of the 100 simulations. The detection of the true causal structures in Models S1 and S2 is mainly disturbed by confounding and mediating variables, respectively, which are typical challenges for causal inference. In Model S1, all methods correctly identified the causal links \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e, and \u003cem\u003eS\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e. The CCF method incorrectly identified the links \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e due to the influence of the confounding variables \u003cem\u003eP\u003c/em\u003e (simultaneously affecting \u003cem\u003eS\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e) and \u003cem\u003eS\u003c/em\u003e (simultaneously affecting \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eI\u003c/em\u003e) respectively. Besides, the indirect link \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e results from the mediating variable \u003cem\u003eS\u003c/em\u003e. These incorrect links in other causal methods (CCM, TE, and PCMCI+) also share the similar reasons. As the bivariate unidirectional causal methods, the CCF and TE methods miss the link \u003cem\u003eI\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e because the opposite direction \u003cem\u003eS\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e has a more significant causal effect. It should be noted that since the CCM method successfully identified all of the predefined causal links, many new false connections, i.e., \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e, \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e, \u003cem\u003eS\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e, \u003cem\u003eI\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e, were generated with the control of confounding and mediating variables. One reasonable explanation is that the synthetic system is strongly coupled, such that the CCM method easily mistakes unidirectional causal relationships for bidirectional relationships. In Model S2, PCMCI\u0026thinsp;+\u0026thinsp;perfectly restores all predefined causal links, while the other three methods (CCF, CCM, and TE) incorrectly identified many indirect causal links due to the transitivity of causal relationships. Similar to the results of Model 1, the CCF and TE methods miss the causal link \u003cem\u003eI\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e, and CCM also mistakes all unidirectional links for bidirectional ones.\u003c/p\u003e\n \u003cp\u003eTable \u003cspan\u003e3\u003c/span\u003e further lists the evaluation results of the causal tests in the synthetic case. The PCMCI\u0026thinsp;+\u0026thinsp;method shows the best overall performance in both model structures, with the TPR, AR, and FPR values of 0.97(1.00), 0.85(0.96), and 0.23(0.05) in Model S1(S2). The CCF and TE methods show relatively lower TPR and AR values, and relatively higher FPR values. In contrast, the CCM method shows the lowest AR (lower than 0.5) and the highest TPR and FPR (higher than 0.9). Besides, in situations where direct and indirect causality are not distinguished, that is, the indirect links identified by the methods are not considered as the wrong links (the values in parentheses of Table \u003cspan\u003e3\u003c/span\u003e), the performances of CCF, CCM, and TE methods improve significantly, especially in Model S2, indicating that the three methods are strongly affected by mediating variables and the identified causal links are likely to contain indirect effects, which might hinder our understanding of the real mechanism in hydrological systems. It should be noted that in such situations, the TPR values are not changed because the true positives (TP, i.e., the causal links generated by the hydrological model are correctly identified) remain the same. In addition, due to the complexity of the hydrological model, the overall performance of all causal methods in Model S2 is better than that in Model S1 except for the CCM method, which shows extremely high TPR and FPR simultaneously owing to the mistakes for bidirectional links. Therefore, in the following real case studies, the identified causal links were revised to unidirectional links, that is, the direction with a significant maximum causal effect was regarded as the only causal direction, thus to avoid misidentification.\u003c/p\u003e\n \u003cdiv\u003e \u0026nbsp;\u003ctable id=\"Tab3\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv\u003eTable 3\u003c/div\u003e\n \u003cdiv\u003e\n \u003cp\u003eEvaluation results of the causal tests in the synthetic case.\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\u0026nbsp;\u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCCF\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCCM\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eTE\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003ePCMCI+\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eTPR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.60\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e1.00\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.60\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.97\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eModel S1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.58 (0.60)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.42 (0.50)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.66 (0.69)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e0.85 (0.87)\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eFPR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.43 (0.40)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e1.00 (1.00)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.30 (0.22)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e0.23 (0.22)\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eTPR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.75\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e1.00\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.75\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e1.00\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eModel S2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eAR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.67 (0.89)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.36 (0.49)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.67 (0.89)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e0.96 (1.00)\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003eFPR\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.38 (0.00)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.95 (0.93)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e0.38 (0.00)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"char\"\u003e\n \u003cp\u003e\u003cstrong\u003e0.05 (0.00)\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003ctfoot\u003e\n \u003ctr\u003e\n \u003ctd colspan=\"6\"\u003eNote: The values in parentheses represent situations where direct and indirect causality are not distinguished, i.e., the indirect links identified by the methods are not considered as the wrong links. The TPR values are not changed in such situations.\u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tfoot\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec10\"\u003e\n \u003ch2\u003e3.3 Performance in sensitivity tests\u003c/h2\u003e\n \u003cp\u003eIn this section, the sensitivity of each method in sample length, noise, and missing data was analyzed. Figure \u003cspan\u003e4\u003c/span\u003e presents the impact of sample length with the sizes of 100, 300, 500, 1000, and 2000. 100 datasets were generated to ensure the robustness of the results. As shown in Figs. \u003cspan\u003e4\u003c/span\u003ea and b, in both Model S1 and S2, for CCF and CCM methods, the TPR is not sensitive to variations in sample length; while for TE and PCMCI\u0026thinsp;+\u0026thinsp;methods, the TPR increases with increasing sample length. For the CCF method, which is based on linear lag correlation, 100 samples are enough to detect all the causal links. This number is also applicable to the CCM method, which is based on the deterministic dynamical system theory. In contrast, TE is based on the probability density framework, which needs to estimate the probability density function from the histogram of the frequency distribution, as well as to implement statistical hypothesis testing for conditional independence, thus requiring sufficient sample length. Similarly, the PCMCI\u0026thinsp;+\u0026thinsp;method also requires sufficient length samples in conditional independence tests to keep iterating and removing initialized redundant connections. To achieve relative stability, the TE and PCMCI\u0026thinsp;+\u0026thinsp;methods need at least 1000 and 500 samples for Model S1 and S2, respectively. This difference results from the complexity of the model structure.\u003c/p\u003e\n \u003cp\u003eAs for the Accuracy Rate (AR), the values in each causal method fluctuate or even decrease with the increasing sample length, especially for TE and CCM methods (Figs. \u003cspan\u003e4\u003c/span\u003ec and d). This is attributed to the limitations of causal methods in dealing with the impacts of confounding and mediating variables, which become more striking with the increase of simple length and lead to a significant increase in false positive rate. Additionally, the dashed lines in Fig. \u003cspan\u003e4\u003c/span\u003e represent the situations where direct and indirect causality are not distinguished. The difference in trends between the solid and dashed lines can be further analyzed to determine whether variations in sample length affect the identification of indirect causal links. The CCF, CCM, and PCMCI\u0026thinsp;+\u0026thinsp;methods show a similar trend between the solid and dashed lines in both model structures, while for the TE method, the variation rate is slightly inconsistent or even the reverse, due to the increasing detection rate of indirect links \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e in Model S1, and \u003cem\u003eT\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e, \u003cem\u003eT\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, \u003cem\u003eM\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e in Model S2, with the increase of sample length. Besides, in the situation where direct and indirect causality is not distinguished, the improvement in the performance of causal methods in Model S2 is more pronounced than that in Model S1, which is consistent with the setup of model structure, i.e., the former is mainly controlled by indirect causality.\u003c/p\u003e\n \u003cp\u003eFigure \u003cspan\u003e5\u003c/span\u003e presents the impact of noise with the SNR(dB) of 2, 3, 5, 10, 20, 30, and 40. The sample length was fixed as 2000, and 100 datasets were generated to ensure the robustness of the results. As shown in Fig. \u003cspan\u003e5\u003c/span\u003ea, in model S1, the noise has little impact on the TPR values for the CCF and CCM methods with the increase of noise level (i.e., the decrease of SNR). In contrast, for the TE method, the TPR increases with the increasing noise level, especially when the SNR (dB) is below 20. This is attributed to the assumption of the nonlinear stochastic system for the TE method, and appropriate noise helps to identify causal relationships. For PCMCI+, the TPR decreases with the decreasing SNR(dB) and shows a slight increase when the SNR (dB) is lower than 10. However, in Model S2, the noise has little impact on the TPR values for all causal methods (Fig. \u003cspan\u003e5\u003c/span\u003eb), owing to the insensitivity of this model structure to input noise. As for the Accuracy Rate (AR), the values in CCF and TE are relatively stable with the variations of the noise level in both model structures (Figs. \u003cspan\u003e5\u003c/span\u003ec and d), while for the PCMCI\u0026thinsp;+\u0026thinsp;method, the AR decreases with the decreasing SNR in Model S1, and for the CCM method, the AR increases as the SNR(dB) drops below 10 in both model structures. Additionally, all causal methods show a similar trend between the solid and dashed lines in both model structures, indicating that the noise does not have significant impacts on the identification of indirect causal links.\u003c/p\u003e\n \u003cp\u003eFigure \u003cspan\u003e6\u003c/span\u003e presents the impact of missing data on the performance of each causal method with the sparsity rate increasing from 10\u0026ndash;60%. Two scenarios of missing data, namely synchronous and asynchronous sparsity, which may occur in real hydrologic data due to equipment errors and defective storage technologies, were constructed with 100 tests. To ensure the computability of all causal methods, the missing values were complemented by the arithmetic mean of observed values of the same variable. As shown in Figs. \u003cspan\u003e6\u003c/span\u003ea and b, with the increase of synchronous missing rates, the TPR values remain relatively stable for the CCF method and begin to decline at a critical missing value for the other three methods, especially for CCM and PCMCI+. In comparison with the synchronous sparsity, the TPR values remain relatively stable for the PCMCI method in the scenario of asynchronous sparsity, while for the CCM method, the TPR values decline significantly over the 30% point of missing rate (Figs. \u003cspan\u003e6\u003c/span\u003ee and f). One reasonable explanation is that the CCM method requires state space reconstruction by continuous time series data, while the missing data can lose numerous effective information for the state space, especially in the asynchronous scenario, thus disturbing the detection of causal relationships.\u003c/p\u003e\n \u003cp\u003eAs for the Accuracy Rate (AR), with the increasing missing rate, the values in the CCF method remain relatively stable, while in the other three methods (CCM, TE, and PCMCI+), the values show fluctuation (Figs. \u003cspan\u003e6\u003c/span\u003ec, d, f and g). It should be noted that for the CCM method, the AR values increase substantially, especially over the critical missing value of 30%, which seems implausible. One reasonable explanation is that as the missing rate increases, effective information of variables in their state space is gradually lost, and causally linked variables in the system may no longer maintain the information signature of each other; therefore the efficiency of cross-mapping decreases. The detection rate of both correct and incorrect causal links declines simultaneously, and the latter is more pronounced. In addition, the CCF and PCMCI\u0026thinsp;+\u0026thinsp;methods show a similar trend between the solid and dashed lines, while for the TE method, the variation rate is slightly inconsistent or even the reverse over the 50% point of missing rate due to the decreasing detection rate of indirect links \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e in Model S1, and \u003cem\u003eT\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e and \u003cem\u003eT\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, \u003cem\u003eM\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e in Model S2, with the increasing missing rate. For the CCM method, the variation rate is extremely inconsistent over the 30% point of the missing rate due to the declining detection rate of all indirect links (\u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eI, I\u003c/em\u003e \u0026rarr; \u003cem\u003eQ, T\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e, \u003cem\u003eT\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, \u003cem\u003eM\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e) simultaneously.\u003c/p\u003e\n \u003cp\u003eAt the end of this section, the results of the sensitivity tests can be summarized as follows: In the impact of sample length, the CCF and CCM methods show relatively stable performance, while the TE and PCMCI\u0026thinsp;+\u0026thinsp;methods require at least 500 samples to achieve relative stability; In the impact of noise, the performance of causal models is affected by the model structure, and the TE and PCMCI\u0026thinsp;+\u0026thinsp;methods perform more unstably for Model S1; In the impact of missing data, all causal methods show relatively stable performance within 30% of the missing rate, and the performance of the CCM method deteriorates rapidly over 30%.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"4. Application to real case study","content":"\u003cp\u003eIn this section, two real study cases are presented to further evaluate the applicability of different causal inference methods in complex hydrological systems.\u003c/p\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003e4.1 Application to real case 1\u003c/h2\u003e \u003cdiv id=\"Sec13\" class=\"Section3\"\u003e \u003ch2\u003e4.1.1 Study area and data\u003c/h2\u003e \u003cp\u003eThe Shale Hills Catchment (SHC; Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e(a)), located in central Pennsylvania, USA, is a V-shaped small (0.08 km\u003csup\u003e2\u003c/sup\u003e) forested headwater catchment with comparatively steep slopes and narrow ridges (Guo et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). The surface elevation ranges from 256 to 310 m, with relatively homogenous land cover and lithology, and relatively heterogeneous soil thickness and organic matter content (Lin, \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2006\u003c/span\u003e). The catchment belongs to a temperate continental climate region and the average air temperature varies from \u0026minus;\u0026thinsp;3\u0026deg;C (January) to 22\u0026deg;C (July). The annual average precipitation is around 980 mm, and the monthly distribution is relatively uniform with a small maximum in summer when the rainfall is usually characterized with high intensity and short duration (Jiang et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). The first-order stream of the catchment converges to the Juniata River, and is usually dry during summer months (Liu et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Besides, given that the snowfall mainly occurs from November to April, this period was set as the snow-cover period in this study, while the other period over the water year (May to October) was divided as the snow-free period, thus to explore different causal mechanisms during wet and dry conditions.\u003c/p\u003e \u003cp\u003eThe hydrological series, including discharge (\u003cem\u003eQ\u003c/em\u003e), precipitation (\u003cem\u003eP\u003c/em\u003e), groundwater level (\u003cem\u003eGW\u003c/em\u003e), soil moisture (\u003cem\u003eSM\u003c/em\u003e), and snow depth (\u003cem\u003eSD\u003c/em\u003e), were obtained from the Critical Zone Observatory Data Site (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.czo.psu.edu/data_time_series.html\u003c/span\u003e\u003cspan address=\"http://www.czo.psu.edu/data_time_series.html\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). The dataset is a reanalysis result from the Flux-PIHM model (Shi et al., \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2013\u003c/span\u003e), and has been subject to strict quality control. The temporal resolution is hourly with the period of 2009/11/1\u0026thinsp;~\u0026thinsp;2010/10/31. The min-max normalized values of the time series in snow-free (May to October) and snow-cover (November to April) periods are shown in Figs.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e(b) and (c), respectively.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section3\"\u003e \u003ch2\u003e4.1.2 Results\u003c/h2\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e presents the causal structures identified by different causal inference methods for the Shale Hills Catchment (SHC) during the snow-free period. Generally, all methods can identify the main causal relationships \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eSM\u003c/em\u003e, \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eGW\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, which represent the basic hydrological processes from precipitation to soil water, groundwater and streamflow, respectively, but show differences in other relationships. For the link between \u003cem\u003eGW\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, The CCF method presents a forward direction, while the CCM and TE methods present a backward direction. Only the PCMCI\u0026thinsp;+\u0026thinsp;method shows a bidirectional relationship, which contributes to the understanding of potential interaction processes between groundwater and streamflow. For the link between \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eGW\u003c/em\u003e, both CCF and CCM methods present a lagged forward link, while the PCMCI\u0026thinsp;+\u0026thinsp;method shows a contemporary backward link. Besides, the CCF and CCM methods present a forward link between \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, which is omitted by the other two methods. In terms of causal strength, all methods support a relatively weak association for \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eGW\u003c/em\u003e, while the CCF and CCM methods tend to favor a stronger causal interaction between \u003cem\u003eGW\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, and the TE and PCMCI\u0026thinsp;+\u0026thinsp;methods tend to favor a stronger link \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eSM.\u003c/em\u003e\u003c/p\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e presents the causal structures identified by different causal inference methods for the SHC during the snow-cover period. The snow depth (SD) variable is added to this hydrological system. Similar to the results during the snow-free period, all methods can detect the main causal relationships \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eGW\u003c/em\u003e, and show differences in other relationships between \u003cem\u003eGW\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, and \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eGW\u003c/em\u003e. However, limitations of causal methods are exposed to this more complex hydrological system. The unreasonable causal links \u003cem\u003eGW\u003c/em\u003e \u0026rarr; \u003cem\u003eSD\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eSD\u003c/em\u003e identified by the CCF and CCM methods might be attributed to their shortcomings in dealing with confounding variables. The causal link \u003cem\u003eSD\u003c/em\u003e \u0026rarr; \u003cem\u003eSM\u003c/em\u003e identified by CCF, CCM, and PCMCI\u0026thinsp;+\u0026thinsp;methods is consistent with the snowmelt process. However, distinguished from the negative causal strength in the PCMCI\u0026thinsp;+\u0026thinsp;method, the CCF method shows an incomprehensible positive causal strength, which might be attributed to the influence of the confounding variable \u003cem\u003eP\u003c/em\u003e (simultaneously affecting \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eSD\u003c/em\u003e). In comparison with the snow-free period, all causal methods in the snow-cover period support a weaker association between \u003cem\u003eP\u003c/em\u003e and \u003cem\u003eSM\u003c/em\u003e and a stronger association between \u003cem\u003eP\u003c/em\u003e and \u003cem\u003eGW\u003c/em\u003e except for the TE method (due to the lack of this causal link) during the snow-cover period.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003e4.2 Application to real case 2\u003c/h2\u003e \u003cdiv id=\"Sec16\" class=\"Section3\"\u003e \u003ch2\u003e4.2.1 Study area and data\u003c/h2\u003e \u003cp\u003eThe Chuosijia River Basin (CRB; Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003e(a)) is situated on the southeastern edge of the Tibetan Plateau, China, with a drainage area of 14813 km\u003csup\u003e2\u003c/sup\u003e and elevation ranging from 2440 to 5403m (Yang et al., \u003cspan citationid=\"CR81\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Impacted by the westerly circulation and the southwest monsoon, the monsoon climate here is remarkable with distinct dry and wet seasons, and the same period of rain and heat. The multi-annual average temperature and precipitation are 8.6 ℃ and 740mm, respectively. The intra-annual distribution of precipitation is uneven. More than 80% of the precipitation is concentrated from June to October and usually in the form of heavy rainfall due to the warm moisture from the southwestern Indian Ocean (Chen and Alexander, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). The multi-year average runoff depth is 392mm, and runoff is primarily formed by precipitation, followed by snowmelt and groundwater. The first-order stream of the basin finally converges to Minjiang River, a main tributary of the upper Yangtze River (Liang et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Besides, the CRB has little human interference and can be deemed as a natural basin, which is suitable for investigating causal interactions in natural hydrological systems. Considering that the snowfall mainly occurs from November to May, this period was set as the snow-cover period in this study and the other period over the water year (June to October) was divided as the snow-free period, thus to explore different causal mechanisms during dry and wet conditions.\u003c/p\u003e \u003cp\u003eThe hydrological series includes precipitation (\u003cem\u003eP\u003c/em\u003e), soil moisture (\u003cem\u003eSM\u003c/em\u003e), evaporation (\u003cem\u003eE\u003c/em\u003e), snow water equivalent (\u003cem\u003eSWE\u003c/em\u003e), and discharge (\u003cem\u003eQ\u003c/em\u003e). The precipitation data was obtained from CN05.1(Wu and Gao, \u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e2013\u003c/span\u003e), a gridded dataset (0.25\u0026deg;\u0026times;0.25\u0026deg;; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://ccrc.iap.ac.cn/\u003c/span\u003e\u003cspan address=\"https://ccrc.iap.ac.cn/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) based on more than 2400 observed meteorological stations. The soil moisture data (0-100 cm) was gained from SMCI1.0 (Li et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), a long-term high-resolution soil moisture dataset (1km\u0026times;1km; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://data.tpdc.ac.cn/home/\u003c/span\u003e\u003cspan address=\"https://data.tpdc.ac.cn/home/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) based on the measurements of 1789 observed stations across China. The evaporation and snow water equivalent data were obtained from GLEAM (0.25\u0026deg;\u0026times;0.25\u0026deg;; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.gleam.eu/\u003c/span\u003e\u003cspan address=\"https://www.gleam.eu/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and ERA5-land (0.1\u0026deg;\u0026times;0.1\u0026deg;; \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cds.climate.copernicus.eu/#!/home/\u003c/span\u003e\u003cspan address=\"https://cds.climate.copernicus.eu/#!/home/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), respectively, which have been widely applied in the Tibetan Plateau recently (Chen et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Li et al., \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Yang et al., \u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). The discharge data was gained from the Hydrological Yearbook of the Yangtze River Basin. All gridded data was aggregated to the whole basin scale, with the daily temporal resolution and the whole period from 2008 to 2018. The min-max normalized values of the time series in snow-free (June to October) and snow-cover (November to May) periods during water year 2015\u0026ndash;2016 are shown in Figs.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003e(b) and (c), respectively.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section3\"\u003e \u003ch2\u003e4.2.2 Results\u003c/h2\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003e presents the causal structures identified by different causal inference methods for the Chuosijia River Basin (CRB) during the snow-free period. In general, all methods can identify the main causal relationships \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, which denote the basic hydrological processes from precipitation to soil water and streamflow, respectively, but show differences in other relationships. For the link between \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e, the CCF and CCM methods present a lagged forward link, while the PCMCI\u0026thinsp;+\u0026thinsp;method shows a contemporary backward link. For the link between \u003cem\u003eE\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e, all methods show a forward direction except for the CCM method. For the links between \u003cem\u003eE\u003c/em\u003e and \u003cem\u003eSM\u003c/em\u003e, and \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eE\u003c/em\u003e, the CCF and PCMCI\u0026thinsp;+\u0026thinsp;methods present forward links, while the CCM method shows backward links, providing a complementary understanding of the interactions among daily evaporation, soil moisture, and streamflow. With respect to the causal strength, all methods support a relatively strong association for \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ.\u003c/em\u003e\u003c/p\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003e shows the causal structures identified by different causal inference methods for the CRB during the snow-cover period. The snow water equivalent (\u003cem\u003eSWE\u003c/em\u003e) variable is added to this hydrological system. Generally, all methods can identify the causal relationships \u003cem\u003eE\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e, which reflect the basic processes of the hydrologic cycle, i.e., from evaporation to precipitation and then to streamflow, but show large differences in other relationships. For the link between \u003cem\u003eP\u003c/em\u003e and \u003cem\u003eSM\u003c/em\u003e, the CCF and PCMCI\u0026thinsp;+\u0026thinsp;methods present a forward link, while the CCM and TE methods show a backward link. Besides, the CCF and CCM methods identify the causal links \u003cem\u003eE\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eE\u003c/em\u003e \u0026rarr; \u003cem\u003eSM\u003c/em\u003e, which are omitted by the other two methods. Besides, except for the TE method, all methods identify the causal links \u003cem\u003eSM\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e and \u003cem\u003eE\u003c/em\u003e \u0026rarr; \u003cem\u003eSWE\u003c/em\u003e, which reflect subsurface flow and snowmelt processes, respectively. Only the CCM method detects the snow cover process (\u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eSWE\u003c/em\u003e). Nevertheless, many unreasonable causal relationships emerge in this five-variable system, for instance, the \u003cem\u003eSWE\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e detected by the CCF and TE methods, and the \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eSWE\u003c/em\u003e detected by the CCM method. Compared with the snow-free period, all causal methods in the snow-cover period support a stronger association between \u003cem\u003eE\u003c/em\u003e and \u003cem\u003eP\u003c/em\u003e and a weaker association between \u003cem\u003eP\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e in the snow-cover period, and the CCF and CCM methods support a stronger association between \u003cem\u003eSM\u003c/em\u003e and \u003cem\u003eQ\u003c/em\u003e. The results indicates the increasing contribution of baseflow to streamflow during the snow-free period, as well as the diminishing effect of precipitation on streamflow.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003e5.1 Comparison of four causal inference methods\u003c/h2\u003e \u003cp\u003eIdentifying causal associations in hydrological systems is challenging. Applying different causal methods to artificially generated and real-world hydrological data has been confirmed to yield some incomprehensible or even contradictory results, which leads to some reasonable doubts about the reliability of the currently popular methods, and highlights the importance of maintaining critical attention to the limitations of individual methods.\u003c/p\u003e \u003cp\u003eThe CCF method, which is rooted in traditional linear lag correlation theory, shows numerous significant connections whether in synthetic or real cases. Since reducing the significance level may control the number of false links (Rinderera et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2018\u003c/span\u003e), yet, significant associations can still be detected as the \u003cem\u003ep\u003c/em\u003e-value decreases to 0.001. This might be attributed to relatively strong coupling and synchronization of the hydrological systems. Besides, the essence of the CCF method lies in statistical dependencies, without excluding indirect and confounding factors (Dean and Dunsmuir, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2016\u003c/span\u003e), which leads to many spurious causal associations, such as \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eS\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea), \u003cem\u003eM\u003c/em\u003e \u0026rarr; \u003cem\u003eI\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee) in the synthetic case, and \u003cem\u003eGW\u003c/em\u003e \u0026rarr; \u003cem\u003eSD\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ea), \u003cem\u003eSWE\u003c/em\u003e \u0026rarr; \u003cem\u003eP\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e12\u003c/span\u003ea) in the real cases 1 and 2, respectively. Nevertheless, the CCF method is simple and efficient, and presents a stable performance in the sensitivity tests, with little variations induced by the sample length, noise, and missing values. Therefore, the CCF method can be used as a preliminary means of inferring causal associations in hydrological systems, and a benchmark for the comparison of other causal inference methods.\u003c/p\u003e \u003cp\u003eFor the CCM method, results in the synthetic case show extremely high false positive rates (almost 1) in both model structures, even if all preset causal relationships were identified simultaneously (Figs.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb, f). The irrational bi-directional causality was also reported by Ombadi et al. (\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), which highlights the difficulty for the CCM method to address strong coupling in hydrological systems. In real cases, the causal directions were revised to that with the significant maximum causal effect, while some unreasonable causal links were still identified, such as \u003cem\u003eGW\u003c/em\u003e \u0026rarr; \u003cem\u003eSD\u003c/em\u003e, \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eSD\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eb), and \u003cem\u003eQ\u003c/em\u003e \u0026rarr; \u003cem\u003eSWE\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e12\u003c/span\u003eb) in the real case 1 and 2, respectively, which further confirms the limitations of CCM in dealing with confounding factors (Delforge et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Nevertheless, the CCM method still helps to understand complex hydrological processes, which requires us to return to the basic principles and assumptions of this method, that is, the nonlinear causality identified by CCM would only make sense under the framework of nonlinear dynamic system (NDS) with its concept fully understood (Zhao et al., \u003cspan citationid=\"CR86\" class=\"CitationRef\"\u003e2023a\u003c/span\u003e). The CCM method might contribute to exploring hidden hydrological causal mechanisms from the perspective of NDS, such as the interactions between groundwater and streamflow (Bonotto et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), and within cryospheric hydrological cycle (Zhao et al., \u003cspan citationid=\"CR87\" class=\"CitationRef\"\u003e2023b\u003c/span\u003e), which may not be directly observed from the geophysical level, and thus in turn helps to develop models based on physical foundations.\u003c/p\u003e \u003cp\u003eFor the TE method, the results exhibit better performance in synthetic large sample tests, compared with the CCF and CCM methods, especially in cases where direct and indirect causality are not distinguished. In real cases, the identified significant causal relationships are relatively few, but most of them can be explained by basic hydrological laws, such as the links \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eSM\u003c/em\u003e (Figs.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003ec, \u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ec) and \u003cem\u003eP\u003c/em\u003e \u0026rarr; \u003cem\u003eQ\u003c/em\u003e (Figs.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003ec, \u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e12\u003c/span\u003ec). Based on information theory, the TE method is not constrained by linear assumptions and partly overcomes the autocorrelation of hydrological time series using conditional mutual information, thus can be well applied in complex hydrological systems even with the threshold effect (Moges et al., \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Tennant et al., \u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). However, similar to the CCF and CCM methods, it fails to fundamentally avoid the interference of confounders and indirect causality, and some unreasonable connections still need to be cautiously explained. Besides, the sensitivity analysis indicates the necessity of sufficient sample length (i.e., at least 500\u0026ndash;1000) for TE analysis, due to the challenge of estimating three-dimensional probability density, which has also been emphasized in other research (Kathpalia et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Ruddell and Kumar, \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2009\u003c/span\u003e)\u003c/p\u003e \u003cp\u003eThe PCMCI\u0026thinsp;+\u0026thinsp;method, which is based on causal network learning, shows optimal performance in synthetic large-sample tests, and satisfactory interpretability in real cases. This method infers causality from the network graph of multivariate time series, utilizing iterative conditional independence tests to largely overcome the interference of autocorrelation and confounding factors, and distinguish direct and indirect causality in hydrological systems (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed, h). Notwithstanding, contemporary causal associations in real cases, partly due to the coarse time resolutions (Runge et al., \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2019a\u003c/span\u003e), need to be explained with caution. In terms of the algorithms for the conditional independence test, partial correlation (ParCorr) was selected in this study, which has also been widely adopted in other research (Karmouche et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Nowack et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Since other testing algorithms, such as conditional mutual information (CMI), seems suitable for discovering nonlinear associations in hydrological systems, yet it is fairly time-consuming (Runge et al., \u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e2019b\u003c/span\u003e), and recent research showed a low recall rate in synthetic datasets, and an unstable performance in real cases (Delforge et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), which highlights the importance of choosing appropriate testing algorithms. Moreover, the graph-based method requires the assumption of causal sufficiency, i.e. the absence of hidden variables. In general, this assumption cannot be directly tested, and researchers can only consider known factors, as many as possible, in the causal network analysis through their prior expert knowledge, which may add potential uncertainty to the results. Yet, some recent work indicates that the assumption could be appropriately relaxed, especially in the real world (Gerhardus, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Runge, \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2018\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eIn summary, we recommend a flexible strategy for selecting causal inference methods based on the purpose of research. For exploratory work, such as investigating possible causal connections, the CCF and CCM methods are suggested, even if some irrational links are inevitably produced. The former detects all possible linear associations in hydrological systems, while the latter can reveal the hidden causal interactions in nonlinear dynamic systems. For confirmatory work, such as selecting significant predictors, the TE method can be selected with relatively low false detections. The PCMCI\u0026thinsp;+\u0026thinsp;method, with optimal integrative performance and remarkable interpretability, can be applied to both works mentioned above, and is suggested to deeply understand the interaction mechanisms among multivariates in complex hydrological systems.\u003c/p\u003e \u003cp\u003eNevertheless, each method improves, more or less, our understanding of hydrological processes, which calls into our multi-perspective comprehension of the findings obtained by these methods under the context of their particular assumptions and constraints. Specifically, the causality identified by the CCF method can be understood as the linear connection, while identified by the CCM, TE, and PCMCI\u0026thinsp;+\u0026thinsp;methods need to be understood from nonlinear dynamics, information theory, and causal network perspectives, respectively. Therefore, it is important to combine the priori expert knowledge with the assumptions and limitations of each method, to comprehensively explain the reasonable or implausible causal associations. Moreover, the application of different causal methods may yield conflicting results, which highlights the importance to state the assumptions when drawing conclusions of causal issues. To this end, hydrological researchers are encouraged to be more explicit in elaborating assumptions that enable more robust conclusions, and in interpreting and evaluating conclusions under alternative sets of assumptions (Runge et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2023\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003e5.2 Direct and indirect causality\u003c/h2\u003e \u003cp\u003eAs the central focus of hydrology, the water cycle contains numerous components, including rainfall, evaporation, soil moisture, runoff, groundwater, snow, etc., which are connected in direct (e.g., rainfall-soil moisture) or indirect (e.g. groundwater-evaporation) way. Causal inference contributes to revealing the functional connectivity of these components in hydrological systems (Delforge et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Rinderera et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). However, the causal methods may erroneously identify indirect causal influences as direct ones due to the transitivity of the causal relationship (Park et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). On the other hand, if two components are not connected directly, a causal influence may also exist, but it must be indirect. It is of great importance to remove the effects of indirect causality and thus to determine the direct causal relationships, as the latter serves as the foundation for modeling, prediction, and control of the system (Leng et al., \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). In comparison with other algorithms, the PCMCI\u0026thinsp;+\u0026thinsp;method can effectively distinguish between direct and indirect causal links (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e; Table\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e3\u003c/span\u003e), with explicit indirect impact mechanisms from the graph of multivariate time series. Nevertheless, discriminating between direct and indirect causality for the real hydrological systems is still challenging, because both relationships are often intertwined. Taking the \u003cem\u003eP\u003c/em\u003e-\u003cem\u003eSM-Q\u003c/em\u003e process for example, the precipitation falling on land can first supplement soil moisture through infiltration and then move to the outlet of the basin through subsurface processes (i.e., the indirect process), or directly transfer through surface processes such as saturation overland flow or infiltration excess runoff (i.e., the direct process) (Kidron, \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). How to quantify such complex causal interactions among three or even more variables in hydrological systems remains to be further considered (Goodwell and Bassiouni, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Weijs et al., \u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e2018\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003e5.3 Progress, limitations, and future perspectives\u003c/h2\u003e \u003cp\u003eIn comparison with the previous research (Delforge et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Ombadi et al., \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), our study focuses on hydrological systems, introduces the latest causal network learning algorithm PCMCI+, and serves the CCF method as the link between causality and correlation, and as the benchmark for the comparison of other methods. Moreover, in the synthetic cases, we expanded the structures of the conceptual hydrological models, enriched the sensitivity tests, and discussed the direct and indirect causality in hydrological systems for the first time. Two real cases deepen our understanding of different causal methods at different spatiotemporal scales, respectively. Our research in both cases systematically investigated the question of when and how to apply these methods, and how to interpret their results, thus contributing to their further popularity in the hydrology community.\u003c/p\u003e \u003cp\u003eNonetheless, this study concentrates more on the methodological issues, the inner mechanisms of hydrological processes across spatio-temporal scales in the study area are beyond the scope of this study, which need further research (Liu et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Wen et al., \u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Hao et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Besides, under different causal frameworks, testing causal mechanisms either unduly relies on assumptions or lacks theoretical examination (Su et al., \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Thus, the discovery and test of causal mechanisms in complicated hydrological systems with unknown causal structures is still challenging. Moreover, the hydrological systems may involve certain latent variables, such as unknown/hidden or unobservable variables, which can introduce some confounding factors and make the detected links spurious. Fortunately, the LPCMCI algorithm, proposed by Runge (\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) recently, may overcome this limitation to some extent and can be further applied in hydrology. Additionally, some new methods, such as the Partial cross mapping (PCM; Leng et al., \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) and the method proposed by Park et al. (\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2023\u003c/span\u003e), which address the issue of indirect causality based on nonlinear dynamics theory and monotonic ODE model, respectively, can be tested and evaluated in the future.\u003c/p\u003e \u003cp\u003eRecently, the questions\u0026ndash;assumptions\u0026ndash;data (QAD) template, proposed by Runge et al. (\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2023\u003c/span\u003e), helps to guide researchers on how to phrase and tackle their issues in the framework of causal inference. Yet, some typical challenges, such as the hidden confounding, non-stationarity, contemporaneous causality, and preprocessing of variables, need extra attention. Fortunately, the era of big data provides new opportunities for research in this field, as it could comprehensively characterize the structural information between variables, verify the spurious causal links in the causal network structure, and infer latent variable structures that are difficult to observe (Li et al., 2023). To this end, causal inference helps to build the bridge between data-driven machine learning and prior expert knowledge, thus promoting the understanding of complex hydrological processes and the robust causal prediction.\u003c/p\u003e \u003c/div\u003e"},{"header":"Conclusion","content":"\u003cp\u003eIn this research, the performances of four popular causal inference methods (CCF, CCM, TE, and PCMCI+) were systematically evaluated in hydrological systems using both synthetic and real-world time series. For the synthetic cases, the PCMCI\u0026thinsp;+\u0026thinsp;method shows the best performance in large sample tests, while the other three methods present relatively poor performances due to the interference of confounding and mediating factors. In sensitivity tests, the CCF method is less affected by sample length, noise, and missing values, while the CCM method is significantly impacted by the missing rate, especially over the critical value of 30%. Besides, the TE and PCMCI\u0026thinsp;+\u0026thinsp;methods need at least 500 samples to achieve relative stability, and are vulnerable to the interference of noise. For the real cases, the PCMCI\u0026thinsp;+\u0026thinsp;method shows the best interpretability, while the other three methods produce many inexplicable causal links. Additionally, some contradictory results among different methods emerge, owing to their different assumptions and limitations.\u003c/p\u003e \u003cp\u003eIn summary, the PCMCI\u0026thinsp;+\u0026thinsp;method serves as a favorable choice for conducting causal inference within hydrological systems, while some strong assumptions, such as the causal sufficiency, need to be considered with caution. Nonetheless, each method improves, more or less, our understanding of hydrological processes, which requires our multi-perspective comprehension of the findings obtained by these methods under the context of their particular assumptions and constraints. A comprehensive application of diverse methods based on the specific issue is encouraged for the robustness of conclusions, with the assumptions stated explicitly in advance. Promisingly, the causal inference methods provide a complementary data-driven avenue to unlock the inner mechanisms of complex hydrological systems, and have broad application prospects in the hydrology community.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe research is supported by the National Nature Science Foundation of China (No. 52379017).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCRediT authorship contribution statement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eHanxu Liang: Methodology, Investigation, Formal analysis, Visualization, Writing-original draft. Wensheng Wang: Supervision, Funding acquisition, Writing - Review \u0026amp; Editing. Bin Chen: Data curation. Li Guo: Writing - Review \u0026amp; Editing. Hu Liu: Writing - Review \u0026amp; Editing. Siyi Yu: Validation. Dan Zhang: Conceptualization, Validation.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eDeclaration of Competing Interest\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll authors agreed to the published version of the manuscript and declare no conflicts of interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAltman, N., and M. Krzywinski (2015), Association, correlation and causation, \u003cem\u003eNat. Methods\u003c/em\u003e \u003cem\u003e12\u003c/em\u003e(10), 899-900. https://doi.org/10.1038/nmeth.3587.\u003c/li\u003e\n\u003cli\u003eAngelov, P. P., E. A. Soares, R. C. Jiang, N. I. Arnold, and P. M. Atkinson (2021), Explainable artificial intelligence: an analytical review, \u003cem\u003eWires. Data. Min. Knowl.\u003c/em\u003e \u003cem\u003e11\u003c/em\u003e(5). https://doi.org/10.1002/widm.1424.\u003c/li\u003e\n\u003cli\u003eApaydin, H., and M. Sibtain (2021), A multivariate streamflow forecasting model by integrating improved complete ensemble empirical mode decomposition with additive noise, sample entropy, Gini index and sequence-to-sequence approaches, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e603\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2021.126831.\u003c/li\u003e\n\u003cli\u003eBennett, A., B. Nijssen, G. X. Ou, M. Clark, and G. Nearing (2019), Quantifying Process Connectivity With Transfer Entropy in Hydrologic Models, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e55\u003c/em\u003e(6), 4613-4629. https://doi.org/10.1029/2018wr024555.\u003c/li\u003e\n\u003cli\u003eBergstr\u0026ouml;m, S., and L. P. Graham (1998), On the scale problem in hydrological modelling, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e211\u003c/em\u003e(1-4), 253-265. https://doi.org/10.1016/s0022-1694(98)00248-0.\u003c/li\u003e\n\u003cli\u003eBl\u0026ouml;schl, G., M. F. P. Bierkens, A. Chambel, C. Cudennec, and G. Destouni (2019), Twenty-three unsolved problems in hydrology (UPH) - a community perspective, \u003cem\u003eHydrol. Sci. J.\u003c/em\u003e \u003cem\u003e64\u003c/em\u003e(10), 1141-1158. https://doi.org/10.1080/02626667.2019.1620507.\u003c/li\u003e\n\u003cli\u003eBonotto, G., T. J. Peterson, K. Fowler, and A. W. Western (2022), Identifying Causal Interactions Between Groundwater and Streamflow Using Convergent Cross-Mapping, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e58\u003c/em\u003e(8). https://doi.org/10.1029/2021wr030231.\u003c/li\u003e\n\u003cli\u003eBudakoti, S., T. Chauhan, R. Murtugudde, S. Karmakar, and S. Ghosh (2021), Feedback From Vegetation to Interannual Variations of Indian Summer Monsoon Rainfall, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e57\u003c/em\u003e(5). https://doi.org/10.1029/2020wr028750.\u003c/li\u003e\n\u003cli\u003eChen, J., F. P. Brissette, and R. Leconte (2010), A daily stochastic weather generator for preserving low-frequency of climate variability, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e388\u003c/em\u003e(3-4), 480-490. https://doi.org/10.1016/j.jhydrol.2010.05.032.\u003c/li\u003e\n\u003cli\u003eChen, F., W. T. Crow, P. J. Starks, and D. N. Moriasi (2011), Improving hydrologic predictions of a catchment model via assimilation of surface soil moisture, \u003cem\u003eAdv. Water Resour.\u003c/em\u003e \u003cem\u003e34\u003c/em\u003e(4), 526-536. https://doi.org/10.1016/j.advwatres.2011.01.011.\u003c/li\u003e\n\u003cli\u003eChen, R., M. X. Yang, X. J. Wang, G. N. Wan, and H. Y. Li (2022), Thermal regime variations of the uppermost soil layer in the central Tibetan Plateau, \u003cem\u003eCatena\u003c/em\u003e \u003cem\u003e213\u003c/em\u003e. https://doi.org/10.1016/j.catena.2022.106224.\u003c/li\u003e\n\u003cli\u003eChen, Y., and D. Alexander (2022), Integrated flood risk assessment of river basins: Application in the Dadu river basin, China, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e613\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2022.128456.\u003c/li\u003e\n\u003cli\u003eCryer, J. D., and K. Chan (2008), Time series analysis with applications in R. New York, NY: Springer.\u003c/li\u003e\n\u003cli\u003eDean, R. T., and W. T. M. Dunsmuir (2016), Dangers and uses of cross-correlation in analyzing time series in perception, performance, movement, and neuroscience: The importance of constructing transfer function autoregressive models, \u003cem\u003eBehav. Res. Methods\u003c/em\u003e \u003cem\u003e48\u003c/em\u003e(2), 783-802. https://doi.org/10.3758/s13428-015-0611-2.\u003c/li\u003e\n\u003cli\u003eDelforge, D., O. de Viron, M. Vanclooster, M. Van Camp, and A. Watlet (2022), Detecting hydrological connectivity using causal inference from time series: synthetic and real karstic case studies, \u003cem\u003eHydrol. Earth Syst. Sci.\u003c/em\u003e \u003cem\u003e26\u003c/em\u003e(8), 2181-2199. https://doi.org/10.5194/hess-26-2181-2022.\u003c/li\u003e\n\u003cli\u003eFaybishenko, B. (2017), Detecting dynamic causal inference in nonlinear two-phase fracture flow, \u003cem\u003eAdv. Water Resour.\u003c/em\u003e \u003cem\u003e106\u003c/em\u003e, 111-120. https://doi.org/https://doi.org/10.1016/j.advwatres.2017.02.011.\u003c/li\u003e\n\u003cli\u003eGao, Y. B., C. Merz, G. Lischeid, and M. Schneider (2018), A review on missing hydrological data processing, \u003cem\u003eEnviron. Earth. Sci.\u003c/em\u003e \u003cem\u003e77\u003c/em\u003e(2). https://doi.org/10.1007/s12665-018-7228-6.\u003c/li\u003e\n\u003cli\u003eGeng, S., F. Devries, and I. Supit (1986), A simple method for generating daily rainfall data, \u003cem\u003eAgric. Meteorol.\u003c/em\u003e \u003cem\u003e36\u003c/em\u003e(4), 363-376. https://doi.org/10.1016/0168-1923(86)90014-6.\u003c/li\u003e\n\u003cli\u003eGerhardus, A. a. R., Jakob (2020), High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems, volume 33, pages 12615\u0026ndash;12625. Curran Associates, Inc.\u003c/li\u003e\n\u003cli\u003eGong, W., D. W. Yang, H. V. Gupta, and G. Nearing (2014), Estimating information entropy for hydrological data: One-dimensional case, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e50\u003c/em\u003e(6), 5003-5018. https://doi.org/10.1002/2014wr015874.\u003c/li\u003e\n\u003cli\u003eGood, S. P., D. Noone, and G. Bowen (2015), Hydrologic connectivity constrains partitioning of global terrestrial water fluxes, \u003cem\u003eScience\u003c/em\u003e \u003cem\u003e349\u003c/em\u003e(6244), 175-177. https://doi.org/10.1126/science.aaa5931.\u003c/li\u003e\n\u003cli\u003eGoodwell, A. E., and M. Bassiouni (2022), Source Relationships and Model Structures Determine Information Flow Paths in Ecohydrologic Models, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e58\u003c/em\u003e(9). https://doi.org/10.1029/2021wr031164.\u003c/li\u003e\n\u003cli\u003eGranger, C. W. (1969), Investigating causal relations by econometric models and cross‐spectral methods, \u003cem\u003eEconometrica: Journal of the Econometric Society, 37(3), 424\u0026ndash;438\u003c/em\u003e. https://doi.org/10.2307/1912791.\u003c/li\u003e\n\u003cli\u003eGuo, L., J. Chen, and H. Lin (2014), Subsurface lateral preferential flow network revealed by time-lapse ground-penetrating radar in a hillslope, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e50\u003c/em\u003e(12), 9127-9147. https://doi.org/10.1002/2013wr014603.\u003c/li\u003e\n\u003cli\u003eHao, Y., F. B. Sun, H. Wang, W. B. Liu, Y. J. Shen, Z. Li, and S. J. Hu (2022), Understanding climate-induced changes of snow hydrological processes in the Kaidu River Basin through the CemaNeige-GR6J model, \u003cem\u003eCatena\u003c/em\u003e \u003cem\u003e212\u003c/em\u003e. https://doi.org/10.1016/j.catena.2022.106082.\u003c/li\u003e\n\u003cli\u003eJiang, Y. J., Y. L. Zhang, B. H. Fan, J. H. Wen, H. Liu, C. R. Mello, J. F. Cui, C. Yuan, and L. Guo (2023), Preferential flow influences the temporal stability of soil moisture in a headwater catchment, \u003cem\u003eGeoderma\u003c/em\u003e \u003cem\u003e437\u003c/em\u003e. https://doi.org/10.1016/j.geoderma.2023.116590.\u003c/li\u003e\n\u003cli\u003eKarmouche, S., E. Galytska, J. Runge, G. A. Meehl, A. S. Phillips, K. Weigel, and V. Eyring (2023), Regime-oriented causal model evaluation of Atlantic-Pacific teleconnections in CMIP6, \u003cem\u003eEarth System Dynamics\u003c/em\u003e \u003cem\u003e14\u003c/em\u003e(2), 309-344. https://doi.org/10.5194/esd-14-309-2023.\u003c/li\u003e\n\u003cli\u003eKathpalia, A., P. Manshour, and M. Palus (2022), Compression complexity with ordinal patterns for robust causal inference in irregularly sampled time series, \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cem\u003e12\u003c/em\u003e(1). https://doi.org/10.1038/s41598-022-18288-4.\u003c/li\u003e\n\u003cli\u003eKhatibi, R., B. Sivakumar, M. A. Ghorbani, O. Kisi, K. Kocak, and D. F. Zadeh (2012), Investigating chaos in river stage and discharge time series, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e414\u003c/em\u003e, 108-117. https://doi.org/10.1016/j.jhydrol.2011.10.026.\u003c/li\u003e\n\u003cli\u003eKidron, G. J. (2021), Comparing overland flow processes between semiarid and humid regions: Does saturation overland flow take place in semiarid regions?, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e593\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2020.125624.\u003c/li\u003e\n\u003cli\u003eKinney, J. B., and G. S. Atwal (2014), Equitability, mutual information, and the maximal information coefficient, \u003cem\u003eProc. Natl. Acad. Sci. U. S. A.\u003c/em\u003e \u003cem\u003e111\u003c/em\u003e(9), 3354-3359. https://doi.org/10.1073/pnas.1309933111.\u003c/li\u003e\n\u003cli\u003eKleidon, A., and M. Renner (2013), Thermodynamic limits of hydrologic cycling within the Earth system: concepts, estimates and implications, \u003cem\u003eHydrol. Earth Syst. Sci.\u003c/em\u003e \u003cem\u003e17\u003c/em\u003e(7), 2873-2892. https://doi.org/10.5194/hess-17-2873-2013.\u003c/li\u003e\n\u003cli\u003eKonapala, G., S. C. Kao, and N. Addor (2020), Exploring Hydrologic Model Process Connectivity at the Continental Scale Through an Information Theory Approach, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e56\u003c/em\u003e(10). https://doi.org/10.1029/2020wr027340.\u003c/li\u003e\n\u003cli\u003eLeng, S. Y., H. F. Ma, J. Kurths, Y. C. Lai, W. Lin, K. Aihara, and L. N. Chen (2020), Partial cross mapping eliminates indirect causal influences, \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cem\u003e11\u003c/em\u003e(1). https://doi.org/10.1038/s41467-020-16238-0.\u003c/li\u003e\n\u003cli\u003eLi, Q. L., G. S. Shi, W. Shangguan, V. Nourani, J. D. Li, L. Li, F. N. Huang, Y. Zhang, C. Y. Wang, D. G. Wang, J. X. Qiu, X. J. Lu, and Y. J. Dai (2022), A 1 km daily soil moisture dataset over China using in situ measurement and machine learning, \u003cem\u003eEarth. Syst. Sci. Data\u003c/em\u003e \u003cem\u003e14\u003c/em\u003e(12), 5267-5286. https://doi.org/10.5194/essd-14-5267-2022.\u003c/li\u003e\n\u003cli\u003eLi, X., M. Feng, Y. H. Ran, Y. Su, F. Liu, C. L. Huang, H. F. Shen, Q. Xiao, J. B. Su, S. W. Yuan, and H. D. Guo (2023a), Big Data in Earth system science and progress towards a digital twin, \u003cem\u003eNat. Rev. Earth Env.\u003c/em\u003e https://doi.org/10.1038/s43017-023-00409-w.\u003c/li\u003e\n\u003cli\u003eLi, Z. W., X. L. Xu, and K. L. Wang (2023b), Effects of distribution patterns of karst landscapes on runoff and sediment yield in karst watersheds, \u003cem\u003eCatena\u003c/em\u003e \u003cem\u003e223\u003c/em\u003e. https://doi.org/10.1016/j.catena.2023.106947.\u003c/li\u003e\n\u003cli\u003eLi, X. Y., D. Long, Z. Y. Han, B. R. Scanlon, Z. L. Sun, P. F. Han, and A. Z. Hou (2019), Evapotranspiration Estimation for Tibetan Plateau Headwaters Using Conjoint Terrestrial and Atmospheric Water Balances and Multisource Remote Sensing, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e55\u003c/em\u003e(11), 8608-8630. https://doi.org/10.1029/2019wr025196.\u003c/li\u003e\n\u003cli\u003eLiang, H. X., D. Zhang, W. S. Wang, S. Y. Yu, and S. Nimai (2023), Evaluating future water security in the upper Yangtze River Basin under a changing environment, \u003cem\u003eSci. Total Environ.\u003c/em\u003e \u003cem\u003e889\u003c/em\u003e. https://doi.org/10.1016/j.scitotenv.2023.164101.\u003c/li\u003e\n\u003cli\u003eLin, H. (2006), Temporal stability of soil moisture spatial pattern and subsurface preferential flow pathways in the shale hills catchment, \u003cem\u003eVadose Zone J.\u003c/em\u003e \u003cem\u003e5\u003c/em\u003e(1), 317-340. https://doi.org/10.2136/vzj2005.0058.\u003c/li\u003e\n\u003cli\u003eLiu, H., Y. Yu, W. Z. Zhao, L. Guo, J. T. Liu, and Q. Y. Yang (2020), Inferring Subsurface Preferential Flow Features From a Wavelet Analysis of Hydrological Signals in the Shale Hills Catchment, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e56\u003c/em\u003e(11). https://doi.org/10.1029/2019wr026668.\u003c/li\u003e\n\u003cli\u003eLuo, M., F. Meng, Y. Wang, C. Sa, Y. Duan, Y. Bao, and T. Liu (2023), Quantitative detection and attribution of soil moisture heterogeneity and variability in the Mongolian Plateau, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e621\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2023.129673.\u003c/li\u003e\n\u003cli\u003eMassei, N., J. P. Dupont, B. J. Mahler, B. Laignel, M. Fournier, D. Valdes, and S. Ogier (2006), Investigating transport properties and turbidity dynamics of a karst aquifer using correlation, spectral, and wavelet analyses, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e329\u003c/em\u003e(1), 244-257. https://doi.org/https://doi.org/10.1016/j.jhydrol.2006.02.021.\u003c/li\u003e\n\u003cli\u003eMoges, E., B. L. Ruddell, L. Zhang, J. M. Driscoll, and L. G. Larsen (2022), Strength and Memory of Precipitation\u0026apos;s Control Over Streamflow Across the Conterminous United States, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e58\u003c/em\u003e(3). https://doi.org/10.1029/2021wr030186.\u003c/li\u003e\n\u003cli\u003eNowack, P., J. Runge, V. Eyring, and J. D. Haigh (2020), Causal networks for climate model evaluation and constrained projections, \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cem\u003e11\u003c/em\u003e(1). https://doi.org/10.1038/s41467-020-15195-y.\u003c/li\u003e\n\u003cli\u003eOmbadi, M., P. Nguyen, S. Sorooshian, and K. l. Hsu (2020), Evaluation of Methods for Causal Discovery in Hydrometeorological Systems, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e56\u003c/em\u003e(7). https://doi.org/10.1029/2020wr027251.\u003c/li\u003e\n\u003cli\u003ePark, S. H., S. Ha, and J. K. Kim (2023), A general model-based causal inference method overcomes the curse of synchrony and indirect effect, \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cem\u003e14\u003c/em\u003e(1). https://doi.org/10.1038/s41467-023-39983-4.\u003c/li\u003e\n\u003cli\u003ePatil, S., and M. Stieglitz (2014), Modelling daily streamflow at ungauged catchments: what information is necessary?, \u003cem\u003eHydrol. Process.\u003c/em\u003e \u003cem\u003e28\u003c/em\u003e(3), 1159-1169. https://doi.org/10.1002/hyp.9660.\u003c/li\u003e\n\u003cli\u003ePearl, J. (2009), Causality: Models, Reasoning, and Inference, Cambridge University Press, Cambridge, 2 edn. https://doi.org/10.1017/CBO9780511803161.\u003c/li\u003e\n\u003cli\u003ePeng, S. L., K. Mihara, X. L. Xu, K. Kuramochi, Y. Toma, and R. Hatano (2024), Modeling hydrological processes under Multi-Model projections of climate change in a cold region of Hokkaido, Japan, \u003cem\u003eCatena\u003c/em\u003e \u003cem\u003e234\u003c/em\u003e. https://doi.org/10.1016/j.catena.2023.107605.\u003c/li\u003e\n\u003cli\u003eReichenbach, H. (1956), The direction of time, \u003cem\u003eUniversity of California Press, Berkeley and Los Angeles\u003c/em\u003e.\u003c/li\u003e\n\u003cli\u003eReichstein, M., G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, and Prabhat (2019), Deep learning and process understanding for data-driven Earth system science, \u003cem\u003eNature\u003c/em\u003e \u003cem\u003e566\u003c/em\u003e(7743), 195-204. https://doi.org/10.1038/s41586-019-0912-1.\u003c/li\u003e\n\u003cli\u003eRichardson, C. W. (1981), Stochastic simulation of daily precipitation, temperature, and solar-radiation, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e17\u003c/em\u003e(1), 182-190. https://doi.org/10.1029/WR017i001p00182.\u003c/li\u003e\n\u003cli\u003eRinderera, M., G. Ali, and L. G. Larsen (2018), Assessing structural, functional and effective hydrologic connectivity with brain neuroscience methods: State-of-the-art and research directions, \u003cem\u003eEarth-Sci. Rev.\u003c/em\u003e \u003cem\u003e178\u003c/em\u003e, 29-47. https://doi.org/10.1016/j.earscirev.2018.01.009.\u003c/li\u003e\n\u003cli\u003eRuddell, B. L., and P. Kumar (2009), Ecohydrologic process networks: 1. Identification, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e45\u003c/em\u003e(3). https://doi.org/10.1029/2008wr007279.\u003c/li\u003e\n\u003cli\u003eRudin, C. (2019), Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead, \u003cem\u003eNat. Mach. Intell.\u003c/em\u003e \u003cem\u003e1\u003c/em\u003e(5), 206-215. https://doi.org/10.1038/s42256-019-0048-x.\u003c/li\u003e\n\u003cli\u003eRunge, J. (2018), Causal network reconstruction from time series: From theoretical assumptions to practical estimation, \u003cem\u003eChaos\u003c/em\u003e \u003cem\u003e28\u003c/em\u003e(7). https://doi.org/10.1063/1.5025050.\u003c/li\u003e\n\u003cli\u003eRunge, J. (2020), Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Proc. 36th Conf. Uncertainty in Artificial Intelligence (UAI) Vol. 124 of Proc. Machine Learning Research (eds Peters, J. \u0026amp; Sontag, D.) 1388\u0026ndash;1397. .\u003c/li\u003e\n\u003cli\u003eRunge, J., S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Munoz-Mari, E. H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Scholkopf, P. Spirtes, G. Sugihara, J. Sun, K. Zhang, and J. Zscheischler (2019a), Inferring causation from time series in Earth system sciences, \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cem\u003e10\u003c/em\u003e. https://doi.org/10.1038/s41467-019-10105-3.\u003c/li\u003e\n\u003cli\u003eRunge, J., A. Gerhardus, G. Varando, V. Eyring, and G. Camps-Valls (2023), Causal inference for time series, \u003cem\u003eNat. Rev. Earth Env.\u003c/em\u003e \u003cem\u003e4\u003c/em\u003e(7), 487-505. https://doi.org/10.1038/s43017-023-00431-y.\u003c/li\u003e\n\u003cli\u003eRunge, J., P. Nowack, M. Kretschmer, S. Flaxman, and D. Sejdinovic (2019b), Detecting and quantifying causal associations in large nonlinear time series datasets, \u003cem\u003eSci. Adv.\u003c/em\u003e \u003cem\u003e5\u003c/em\u003e(11). https://doi.org/10.1126/sciadv.aau4996.\u003c/li\u003e\n\u003cli\u003eSang, Y. F., V. P. Singh, J. Wen, and C. M. Liu (2015), Gradation of complexity and predictability of hydrological processes, \u003cem\u003eJ. Geophys. Res. Atmos.\u003c/em\u003e \u003cem\u003e120\u003c/em\u003e(11), 5334-5343. https://doi.org/10.1002/2014jd022844.\u003c/li\u003e\n\u003cli\u003eSchreiber, T. (2000), Measuring information transfer, \u003cem\u003ePhys. Rev. Lett.\u003c/em\u003e \u003cem\u003e85\u003c/em\u003e(2), 461-464. https://doi.org/10.1103/PhysRevLett.85.461.\u003c/li\u003e\n\u003cli\u003eShannon, C. E. (1948), A mathematical theory of communication, \u003cem\u003eBell System Technical Journal, 27(3), 379\u0026ndash;423\u003c/em\u003e. https://doi.org/10.1002/j.1538-7305.1948.tb01338.x.\u003c/li\u003e\n\u003cli\u003eShi, Y., K. J. Davis, C. J. Duffy, and X. Yu (2013), Development of a Coupled Land Surface Hydrologic Model and Evaluation at a Critical Zone Observatory, \u003cem\u003eJ. Hydrometeorol.\u003c/em\u003e \u003cem\u003e14\u003c/em\u003e(5), 1401-1420. https://doi.org/10.1175/jhm-d-12-0145.1.\u003c/li\u003e\n\u003cli\u003eShojaie, A., and E. B. Fox (2022), Granger Causality: A Review and Recent Advances, \u003cem\u003eAnnu. Rev. Stat. Appl.\u003c/em\u003e \u003cem\u003e9\u003c/em\u003e, 289-319. https://doi.org/10.1146/annurev-statistics-040120-010930.\u003c/li\u003e\n\u003cli\u003eSivakumar, B. (2004), Chaos theory in geophysics: past, present and future, \u003cem\u003eChaos Solitons \u0026amp; Fractals\u003c/em\u003e \u003cem\u003e19\u003c/em\u003e(2), 441-462. https://doi.org/10.1016/s0960-0779(03)00055-9.\u003c/li\u003e\n\u003cli\u003eStudent (1908), The Probable Error of a Mean, \u003cem\u003eBiometrika\u003c/em\u003e \u003cem\u003e6\u003c/em\u003e(1), 1-25. https://doi.org/10.2307/2331554.\u003c/li\u003e\n\u003cli\u003eSu, J. B., D. X. Chen, D. H. Zheng, Y. Su, and X. Li (2023), The insight of why: Causal inference in Earth system science, \u003cem\u003eScience China-Earth Sciences\u003c/em\u003e. https://doi.org/10.1007/s11430-023-1148-7.\u003c/li\u003e\n\u003cli\u003eSugihara, G., E. R. Deyle, and H. Ye (2017), Misconceptions about causation with synchronyand seasonal drivers reply, \u003cem\u003eProc. Natl. Acad. Sci. U. S. A.\u003c/em\u003e \u003cem\u003e114\u003c/em\u003e(12), E2272-E2274. https://doi.org/10.1073/pnas.1700998114.\u003c/li\u003e\n\u003cli\u003eSugihara, G., R. May, H. Ye, C. H. Hsieh, E. Deyle, M. Fogarty, and S. Munch (2012), Detecting Causality in Complex Ecosystems, \u003cem\u003eScience\u003c/em\u003e \u003cem\u003e338\u003c/em\u003e(6106), 496-500. https://doi.org/10.1126/science.1227079.\u003c/li\u003e\n\u003cli\u003eTakens, F. (1981), Detecting Strange Attractors in Turbulence, \u003cem\u003eIn Dynamical systems and turbulence, Warwick 1980, (pp. 366\u0026ndash;381). Berlin, Heidelberg: Springer\u003c/em\u003e. https://doi.org/10.1007/BFb0091924.\u003c/li\u003e\n\u003cli\u003eTennant, C., L. Larsen, D. Bellugi, E. Moges, L. Zhang, and H. X. Ma (2020), The Utility of Information Flow in Formulating Discharge Forecast Models: A Case Study From an Arid Snow-Dominated Catchment, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e56\u003c/em\u003e(8). https://doi.org/10.1029/2019wr024908.\u003c/li\u003e\n\u003cli\u003eWang, Q. J., C. F. Yue, X. Q. Li, P. Liao, and X. Y. Li (2023), Enhancing robustness of monthly streamflow forecasting model using embedded-feature selection algorithm based on improved gray wolf optimizer, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e617\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2022.128995.\u003c/li\u003e\n\u003cli\u003eWang, Y., J. Yang, Y. Chen, P. De Maeyer, Z. Li, and W. Duan (2018), Detecting the Causal Effect of Soil Moisture on Precipitation Using Convergent Cross Mapping, \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cem\u003e8\u003c/em\u003e. https://doi.org/10.1038/s41598-018-30669-2.\u003c/li\u003e\n\u003cli\u003eWeijs, S. V., H. Foroozand, and A. Kumar (2018), Dependency and Redundancy: How Information Theory Untangles Three Variable Interactions in Environmental Data, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e54\u003c/em\u003e(10), 7143-7148. https://doi.org/10.1029/2018wr022649.\u003c/li\u003e\n\u003cli\u003eWen, H., S. L. Brantley, K. J. Davis, J. M. Duncan, and L. Li (2021), The Limits of Homogenization: What Hydrological Dynamics can a Simple Model Represent at the Catchment Scale?, \u003cem\u003eWater Resour. Res.\u003c/em\u003e \u003cem\u003e57\u003c/em\u003e(6). https://doi.org/10.1029/2020wr029528.\u003c/li\u003e\n\u003cli\u003eWu, J., and X. J. Gao (2013), A gridded daily observation dataset over China region and comparison with the other datasets, \u003cem\u003eChinese Journal of Geophysics-Chinese Edition\u003c/em\u003e \u003cem\u003e56\u003c/em\u003e(4), 1102-1111. https://doi.org/10.6038/cjg20130406.\u003c/li\u003e\n\u003cli\u003eXu, Z. P., X. L. Man, L. L. Duan, and T. J. Cai (2022), Improved subsurface soil moisture prediction from surface soil moisture through the integration of the (de)coupling effect, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e608\u003c/em\u003e, 12. https://doi.org/10.1016/j.jhydrol.2022.127634.\u003c/li\u003e\n\u003cli\u003eYang, W. J., Y. B. Wang, X. Liu, H. P. Zhao, R. Shao, and G. X. Wang (2020), Evaluation of the rescaled complementary principle in the estimation of evaporation on the Tibetan Plateau, \u003cem\u003eSci. Total Environ.\u003c/em\u003e \u003cem\u003e699\u003c/em\u003e. https://doi.org/10.1016/j.scitotenv.2019.134367.\u003c/li\u003e\n\u003cli\u003eYang, Y., S. J. Chen, Y. R. Zhou, G. W. Ma, W. B. Huang, and Y. M. Zhu (2023), Method for quantitatively assessing the impact of an inter-basin water transfer project on ecological environment-power generation in a water supply region, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e618\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2023.129250.\u003c/li\u003e\n\u003cli\u003eYe, H., E. R. Deyle, L. J. Gilarranz, and G. Sugihara (2015), Distinguishing time-delayed causal interactions using convergent cross mapping, \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cem\u003e5\u003c/em\u003e. https://doi.org/10.1038/srep14750.\u003c/li\u003e\n\u003cli\u003eYu, C., B. Gao, R. Mu\u0026ntilde;oz-Carpena, Y. Tian, L. Wu, and O. Perez-Ovilla (2011), A laboratory study of colloid and solute transport in surface runoff on saturated soil, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e402\u003c/em\u003e(1), 159-164. https://doi.org/https://doi.org/10.1016/j.jhydrol.2011.03.011.\u003c/li\u003e\n\u003cli\u003eZhang, D., W. S. Wang, S. Y. Yu, S. Q. Liang, and Q. F. Hu (2021a), Assessment of the Contributions of Climate Change and Human Activities to Runoff Variation: Case Study in Four Subregions of the Jinsha River Basin, China, \u003cem\u003eJ. Hydro. Eng.\u003c/em\u003e \u003cem\u003e26\u003c/em\u003e(9). https://doi.org/10.1061/(asce)he.1943-5584.0002119.\u003c/li\u003e\n\u003cli\u003eZhang, L., E. Moges, J. W. Kirchner, E. Coda, T. C. Liu, A. S. Wymore, Z. X. Xu, and L. G. Larsen (2021b), CHOSEN: A synthesis of hydrometeorological data from intensively monitored catchments and comparative analysis of hydrologic extremes, \u003cem\u003eHydrol. Process.\u003c/em\u003e \u003cem\u003e35\u003c/em\u003e(11). https://doi.org/10.1002/hyp.14429.\u003c/li\u003e\n\u003cli\u003eZhao, Y. Y., T. J. Zhu, Z. Q. Zhou, H. J. Cai, and Z. D. Cao (2023a), Detecting nonlinear information about drought propagation time and rate with nonlinear dynamic system and chaos theory, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e623\u003c/em\u003e. https://doi.org/10.1016/j.jhydrol.2023.129810.\u003c/li\u003e\n\u003cli\u003eZhao, Y. Y., Y. G. Zou, E. Z. Ma, Z. Q. Zhou, Y. Q. Feng, Z. D. Cao, H. J. Cai, C. Li, and Y. H. Yan (2023b), Can groundwater storage in turn affect the cryospheric variables? A new perspective from nonlinear dynamic causality detection, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e624\u003c/em\u003e, 14. https://doi.org/10.1016/j.jhydrol.2023.129910.\u003c/li\u003e\n\u003cli\u003eZounemat-Kermani, M., O. Batelaan, M. Fadaee, and R. Hinkelmann (2021), Ensemble machine learning paradigms in hydrology: A review, \u003cem\u003eJ. Hydrol.\u003c/em\u003e \u003cem\u003e598\u003c/em\u003e, 13. https://doi.org/10.1016/j.jhydrol.2021.126266.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"stochastic-environmental-research-and-risk-assessment","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"serr","sideBox":"Learn more about [Stochastic Environmental Research and Risk Assessment](https://www.springer.com/journal/477)","snPcode":"477","submissionUrl":"https://submission.nature.com/new-submission/477/3","title":"Stochastic Environmental Research and Risk Assessment","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"causal inference, cross-correlation function (CCF), convergent cross mapping (CCM), transfer entropy (TE), PCMCI+, hydrological systems","lastPublishedDoi":"10.21203/rs.3.rs-4643196/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4643196/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eMany research issues in hydrological systems are intrinsically causal, aiming to determine whether and how one factor affects another. Although causal inference methods have been applied more or less in hydrology, there still remains a lack of systematic comparison between different methods. Here, four popular methods in the causal inference community, including the cross-correlation function (CCF), convergent cross mapping (CCM), transfer entropy (TE), and a causal network learning algorithm (PCMCI+) were selected, with a detailed explanation of their basic principles and underlying assumptions. Next, the performances of these methods were evaluated in large sample tests and sensitivity analysis using synthetic time series generated by a conceptual hydrological model with two predesigned causal structures. Then, the four methods were applied in two real-world cases to further understand their characteristics. The findings show the superior performance of the PCMCI\u0026thinsp;+\u0026thinsp;method in synthetic cases and a commendable level of interpretability in real cases, thus warranting its broader application in hydrological systems. The limitations of the other three methods, especially in effectively addressing confounding and mediating factors, led to several unreasonable causal links. Furthermore, the emergence of conflicting results among different methods in real-world applications underscores the necessity for a multifaceted understanding based on their particular assumptions and constraints. A comprehensive application of diverse methods according to the specific issue is encouraged for the robustness of conclusions, with their assumptions clearly stated in advance. Overall, our research reveals the potential and limitations of different causal inference methods in comprehension of complex interactions within hydrological systems, serving as a useful guide for their further prosperity in hydrology.\u003c/p\u003e","manuscriptTitle":"Inferring causal associations in hydrological systems: A comparison of methods","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-07-19 19:28:01","doi":"10.21203/rs.3.rs-4643196/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2024-09-23T17:08:28+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2024-08-26T13:51:03+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2024-08-14T20:00:38+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"117723746444702419050862116238714957434","date":"2024-07-30T06:24:18+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"160989699627236138633854131429955206559","date":"2024-07-08T14:56:04+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2024-07-07T23:35:36+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2024-06-28T10:44:15+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2024-06-26T14:48:35+00:00","index":"","fulltext":""},{"type":"submitted","content":"Stochastic Environmental Research and Risk Assessment","date":"2024-06-26T13:51:11+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"stochastic-environmental-research-and-risk-assessment","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"serr","sideBox":"Learn more about [Stochastic Environmental Research and Risk Assessment](https://www.springer.com/journal/477)","snPcode":"477","submissionUrl":"https://submission.nature.com/new-submission/477/3","title":"Stochastic Environmental Research and Risk Assessment","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"9157b297-3568-44eb-9ae0-b132e6c5dc08","owner":[],"postedDate":"July 19th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2025-04-21T16:00:13+00:00","versionOfRecord":{"articleIdentity":"rs-4643196","link":"https://doi.org/10.1007/s00477-025-02977-3","journal":{"identity":"stochastic-environmental-research-and-risk-assessment","isVorOnly":false,"title":"Stochastic Environmental Research and Risk Assessment"},"publishedOn":"2025-04-16 15:57:25","publishedOnDateReadable":"April 16th, 2025"},"versionCreatedAt":"2024-07-19 19:28:01","video":"","vorDoi":"10.1007/s00477-025-02977-3","vorDoiUrl":"https://doi.org/10.1007/s00477-025-02977-3","workflowStages":[]},"version":"v1","identity":"rs-4643196","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4643196","identity":"rs-4643196","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}