Causal Reasoning : Towards Dynamic Predictive Models for Runo ff Temporal Behavior of High Dependence Rivers

Nowadays, a noteworthy temporal alteration of traditional hydrological patterns is being observed, producing a higher variability and more unpredictable extreme events worldwide. This is largely due to global warming, which is generating a growing uncertainty over water system behavior, especially river runoff. Understanding these modifications is a crucial and not trivial challenge that requires new analytical strategies like Causality, addressed by Causal Reasoning. Through Causality over runoff series, the hydrological memory and its logical time-dependency structure have been dynamically/stochastically discovered and characterized. This is done in terms of the runoff dependence strength over time. This has allowed determining and quantifying two opposite temporal-fractions within runoff: Temporally Conditioned/Non-conditioned Runoff (TCR/TNCR). Finally, a successful predictive model is proposed and applied to an unregulated stretch, Mijares river catchment (Jucar river basin, Spain), with a very high time-dependency behavior. This research may have important implications over the knowledge of historical rivers ́ behavior and their adaptation. Furthermore, it lays the foundations for reaching an optimum reservoir dimensioning through the building of predictive models of runoff behavior. Regarding reservoir capacity, this research would imply substantial economic/environmental savings. Also, a more sustainable management of river basins through more reliable control reservoirs’ operation is expected to be achieved.

Not all the reasons that explain this increasing variability are brand new.Because of that, there is a consequently strong need to have powerful and reliable analytical methods to build accurate models that reproduce and forecast the future hydrological behavior of a river system [19][20][21][22][23][24][25].Also, there is a growing necessity to design analytical strategies that allow: (a) an increase of knowledge on temporal Water 2019, 11, 877 3 of 18 Indeed, the annual average rainfall in Spain presents a latitudinal decrease pattern, from wet north-west (around 2000 mm), to dry south-east with less than 200 mm [28,38].This is aggravated by low capacity of retention of Spanish soils; on average a 9% & 34% in the rest of Europe [39].
The selected case study was Mijares river (belonging to Jucar river basin, the second largest basin on the Mediterranean side of the Iberian Peninsula).This river is situated in Eastern Spain, traditionally known as "dry Spain".This zone regularly suffers noteworthy drought conditions [40].In particular, an unregulated stretch, characterized by a very high time-dependency behavior and by the existence of important water springs at its headwater was selected [41] (Figure 1a).This case study was defined by gauging station "El Terde", code number 8030 (Figure 1b).
Morphometrically, the sub-basin covers a medium-sized extension of 665 km 2 [42], with a main stream of 47.462 km in length, an average slope of 0.018 meter/meter.Therefore, it has a time of concentration of 12.13 h, according to Spanish regulations [43].
In lithological terms, there are important formations of loams and clays with alternation of gypsum and conglomerates or limestones and gypsums, being that these areas are generally of low permeability.However, these areas may contain deeper aquifers with greater permeability and productivity.In addition, there are small areas of limestones and dolomites where very permeable aquifers are located, generally large and productive [44].
Annual runoff time series (Figure 1c) were supplied through the network of gauging stations belonging to the Jucar River Basin Authority [45].Regarding the historical runoff, it is well known in Hydrology, at least in Spain (South Europe), that the "80 effect" which comprises a drastic reduction in Rainfall and Runoff from the year 80.In this particular basin, this effect is clearly shown in that Figure , where there is also an isolated rise for the period 1988-1991.Most of the authors and experts impute this consequence to Climate Change.
Water 2019, 11, x FOR PEER REVIEW 3 of 19 [18,37].Indeed, the annual average rainfall in Spain presents a latitudinal decrease pattern, from wet north-west (around 2000 mm), to dry south-east with less than 200 mm [28,38].This is aggravated by low capacity of retention of Spanish soils; on average a 9% & 34% in the rest of Europe [39].
The selected case study was Mijares river (belonging to Jucar river basin, the second largest basin on the Mediterranean side of the Iberian Peninsula).This river is situated in Eastern Spain, traditionally known as "dry Spain".This zone regularly suffers noteworthy drought conditions [40].In particular, an unregulated stretch, characterized by a very high time-dependency behavior and by the existence of important water springs at its headwater was selected [41] (Figure 1a).This case study was defined by gauging station "El Terde", code number 8030 (Figure 1b).
Morphometrically, the sub-basin covers a medium-sized extension of 665 km 2 [42], with a main stream of 47.462 km in length, an average slope of 0.018 meter/meter.Therefore, it has a time of concentration of 12.13 h, according to Spanish regulations [43].
In lithological terms, there are important formations of loams and clays with alternation of gypsum and conglomerates or limestones and gypsums, being that these areas are generally of low permeability.However, these areas may contain deeper aquifers with greater permeability and productivity.In addition, there are small areas of limestones and dolomites where very permeable aquifers are located, generally large and productive [44].
Annual runoff time series (Figure 1c) were supplied through the network of gauging stations belonging to the Jucar River Basin Authority [45].Regarding the historical runoff, it is well known in Hydrology, at least in Spain (South Europe), that the "80 effect" which comprises a drastic reduction in Rainfall and Runoff from the year 80.In this particular basin, this effect is clearly shown in that Figure , where there is also an isolated rise for the period 1988-1991.Most of the authors and experts impute this consequence to Climate Change.

Methodololy
This research was articulated in four sequential stages.Firstly, in order to characterize the basin behavior, the memory and temporal behavior indicators of the historical time series were obtained, as well as its main statistical parameters.After that, equiprobable synthetic runoff series were generated through a parsimonious and unconditioned ARMA (1,1) model (Stage-1).Next, synthetic series were applied to populate the Causal Reasoning, supported by Bayesian modelling.This stage was crucial because it comprises the discovery and characterization of the logical and non-trivial time-dependency structure that inherently underlies the hydrological series (Stage-2).Then in Stage-3, based on dependence strength over the time, different dynamic management scenarios (time-step by time-step) were generated.Each temporal dependence scenario contains the dynamic and stochastic identification and quantification of the runoff temporal-fractions, TCR and TNCR.Afterward, an in-depth analysis of each of them was done.Finally, in Stage-4, a preliminary implementation of a dynamic predictive runoff model was performed and its reliability was assessed from a probabilistic perspective (Figure 2).Note: Time series are referred according to water year in Spain where the first month is October and the last is September of following year.

Methodololy
This research was articulated in four sequential stages.Firstly, in order to characterize the basin behavior, the memory and temporal behavior indicators of the historical time series were obtained, as well as its main statistical parameters.After that, equiprobable synthetic runoff series were generated through a parsimonious and unconditioned ARMA (1,1) model (Stage-1).Next, synthetic series were applied to populate the Causal Reasoning, supported by Bayesian modelling.This stage was crucial because it comprises the discovery and characterization of the logical and non-trivial time-dependency structure that inherently underlies the hydrological series (Stage-2).Then in Stage-3, based on dependence strength over the time, different dynamic management scenarios (time-step by time-step) were generated.Each temporal dependence scenario contains the dynamic and stochastic identification and quantification of the runoff temporal-fractions, TCR and TNCR.Afterward, an in-depth analysis of each of them was done.Finally, in Stage-4, a preliminary implementation of a dynamic predictive runoff model was performed and its reliability was assessed from a probabilistic perspective (Figure 2).

Stage-1. Historical Series. Memory and Temporal Behavior Indicators
In order to preliminarily assess the temporal nature of the basin, the memory and temporal behavior indicators of the historical time series were calculated.This was done through Hurst coefficient (H) and correlogram, respectively, based on the general framework exposed in Salas et al. [46] (p.41-44 and p. 49 respectively), and which are well established in the scientific community.Given that the determination of H is not a matter of this research, the reader is referred to Tyralis and Koutsoyiannis [47] where a full discussion about its accuracy is done through three different approaches.After that, main statistical parameters of the historical time series were determined.This essentially comprises getting the basic statistics parameters of the ARMA model, such as mean, standard deviation, and variation and skewness coefficients.
Secondly, a set of equiprobable synthetic runoffs were generated through a parsimonious and unconditioned ARMA (1,1) according to Molina et al. [28] and Molina and Zazo [18,29], based on the classical framework to generate ARMA models as shown by Salas et al. [46] (Chapter 5).A 20 years

Stage-1. Historical Series. Memory and Temporal Behavior Indicators
In order to preliminarily assess the temporal nature of the basin, the memory and temporal behavior indicators of the historical time series were calculated.This was done through Hurst coefficient (H) and correlogram, respectively, based on the general framework exposed in Salas et al. [46] (pp.41-44 and p. 49 respectively), and which are well established in the scientific community.Given that the determination of H is not a matter of this research, the reader is referred to Tyralis and Koutsoyiannis [47] where a full discussion about its accuracy is done through three different approaches.After that, main statistical parameters of the historical time series were determined.This essentially comprises getting the basic statistics parameters of the ARMA model, such as mean, standard deviation, and variation and skewness coefficients.
Secondly, a set of equiprobable synthetic runoffs were generated through a parsimonious and unconditioned ARMA (1,1) according to Molina et al. [28] and Molina and Zazo [18,29], based on the classical framework to generate ARMA models as shown by Salas et al. [46] (Chapter 5).A 20 years warm-up period was defined for achieving a "non-conditioned" synthetic series generation process without "boundary effects".
Water 2019, 11, 877 5 of 18 In the context of this work, it is worth mentioning that an ARMA (p, q) model is formed by adding a moving average (MA) component to the Autoregressive (AR) component to model the temporal correlation in time series [46].AR component ("p" part) represents the temporal dependence delay within a series; the MA component uses delays of the forecast errors to improve the process [28].In this sense, a generic ARMA model is expressed as [48]: where Xt is the value of the variable at a certain time-step t, p is the number of autoregressive parameters, q is the number of moving average parameters, ∅ j and θ j are the coefficient of autoregressive and moving average model respectively and a t is a random variable that represents the historical residuals (error term).

Stage-2. Causal Reasoning
Causality was addressed from the perspective of Causal Reasoning, which is characterized by searching reasoning patterns focused on the cause, and where the objective comprises the prediction of the effect [49].As an analytical strategy, Bayesian modelling was selected, supported by software HUGIN Expert version 7. 3 [50].It should be noted that Bayesian modelling is a powerful AI technique based on a Probabilistic Graphical Model (PGM; [51][52][53][54]) that has significant advantages such as no need for priori information of the process and use of raw data [55,56], its usefulness for analyzing non-linear physical systems [57], or the ease of defining relationships in complex systems and offering a compact representation of the joint probability distribution over sets of random variables [21].
Causal Reasoning is carried out over a set of random decision variables called "nodes" (water years in this case), which are consecutively interconnected through "links" and a set of conditional probability tables between decision variables [21,53,54].Furthermore, the quantification of the variables' relationship strength is performed by Bayes' theorem, which is propagated over time by the conditional probability; in this way the probability distributions are calculated for each decision variable as [58]: where P(A|B) is a probability on event A, assuming that event B is true; P(A,B) is the joint probability of events A and B; and P(B) is the probability of B.
The propagation of the Bayes' Theorem allows, on one hand, dynamically analyze (time-step by time-step) the mitigation (evolution) of temporal dependence over time, through the relative percentage of runoff change that a time-step dynamically produces on the following ones.On the other hand, it allows for generating the yearly interactions among one target year and the previous ones through graphing of marginal dependencies in which each detected connection represents a time-dependency relationship [29].
Hydrological memory of the basin was characterized in terms of the runoff dependence strength over time.This was dynamically performed through interactions over the target year according to the time-horizon defined by the propagation of temporal dependence.The strength of each yearly interaction was evaluated by Dependence Rate (DR; from 0 to 1), expressed as [29]: where p-value is the measure of the strength of evidence against the independence relationship [28,30].Therefore, a DR value equal to 1 expresses a lot of evidence of a total dependence, while in contrast 0 represents little evidence against the independence.DR values were arranged in matrix form and classified in two classes and grouped in 5 intervals (Table 1).This approach, based on a parsimonious and unconditioned ARMA (1,1) combined with the analytical power of Bayesian modeling, confers a high degree of freedom to the Causal Reasoning that overcoming the temporal imposition of the ARMA model [30].This makes it possible to discover non-trivial relationships (time lag > 1) among decision variables, initially consecutively connected.In addition it allows us to discover the hidden logical time-dependency structure, which inherently underlies hydrological series [18,28,29].
For a complete definition of these processes, the reader is invited to refer to Molina et al. [28] and Molina and Zazo [29].

Stage-3. Temporal Runoff Fractions
Once the dependence relationships among variables and propagation of the temporal dependence were discovered through Causality, different management scenarios were dynamically/stochastically generated based on methodological framework shown in Molina and Zazo [18], and a detailed analysis of the temporal behavior was done.
In each time-step (time-lag) of the dependence propagation, both temporal-fractions, TCR (fraction due to time) and TNCR (uninfluenced by time) were determined and quantified.In this sense, the sum of both of them represents the total runoff for each water year.Both fractions are expressed as: where TCR i is the weighted average conditioned runoff of the Target Year (i); w j expresses in weight form the dependence of year j on Target Year (i; please see Figure 7b); Y j is the conditioned runoff for the year j; TD is the temporal interval, expressed in years, of the propagation of dependence from lag = 0 ([0,1] when the relative percentage change is maximum) to maximum lag (when the relative percentage change is minimum and equal to zero; please see Figure 5); R i is the runoff of the time series in the target year; and TNCR i is the uninfluenced by time fraction for the target year.Furthermore, in order to highlight the analytical capacity and discovery potential that this methodological framework offers for the basin water resources management, for each time-lag, two different time series were obtained, one for each temporal fraction.Next, the main statistical parameters and the coefficients of different ARMA models for each time series were determined.Consequently, the differences in the parameters and the coefficients will be exclusively due to "discovered nuances" into the temporal behavior.This was done supported by software TRASERO version 2.3.0 [48].Akaike Information Criterion (AIC; [59]) was applied as criterion of models' selection.This was done by maximizing the expected log-likelihood of a particular model through the maximum likelihood method.The optimal model will be the one with the lowest AIC value [60].This represents the last stage of the research.It comprises, on one side, a preliminary design of the predictive model from a probabilistic perspective, and, on the other, the assessment of the prediction reliability.The predictive model was based on the hidden temporal conditionality that inherently underlies of the historical series and that has been discovered through Causality (Stage-3).
Additionally, some studies have demonstrated the validity of stationary assumptions and stochastic parametric models applied to the generation of hydrological predictive models, such as references [61,62].Recently, Papacharalampous et al. [61] presented an exhaustive study in which large-scale results versus traditional approaches focused on cases studies are provided.Its main conclusion is that stochastic and machine learning methods provide similar forecasts.Earlier, Tyralis and Koutsoyiannis [62] justified the validity of the stationarity hypothesis applied to stochastic methods in a context of prediction of hydroclimatic variables.
In agreement with the above conclusions, the predictive model was performed through stochastic hydrological parametric ARMA models.This was generated by parsimonious ARMA models, one for each temporal fraction (TCR and TNCR).They were developed from a time series obtained in Stage-3, according to temporal horizons of dependence mitigation.
Given that the availability of water resources in Spain is characterized by a high seasonal and annual variability [63][64][65], the predictive model was generated by a time-horizon of one ( 1) year (Dep-1).For that reason, the training dataset was 69 years (time period: 1945/46 to 2013/14) and the prediction was done over water year 2014/15.The obtained value of this final year was used to assess the predictive reliability based on the resulting probability.
Finally, the reliability of the predictive model was assessed from a probabilistic perspective, by suggesting two alternatives.The first one general and the second one detailed.The second approach is based on the fact that the TNCR fraction does not depend on time, and, consequently, adds uncertainty to the prediction.For that, two Monte Carlo simulations from ARMA models with 5000 data were done.The probabilistic approach was based on Gringorten probability expressed as [66]: where F (yi) is the value is the value of the data distribution function y i , i the occupied position by the data in the series ordered from lowest to highest, and n the total number of data in the series.

Statistical Analysis
Figure 3 shows both the main results drawn from the traditional statistical analysis performed (Figure 3a) and the temporal correlation through correlogram (Figure 3b).Please note that the Anderson limits of a correlogram define the temporal non-dependent region; in this way, if the correlation coefficient for a certain time-lag is located inside this region, this time-step is temporally non-dependent.In addition, regarding Hurst coefficient, the resulting value was 0.84, which indicates that analyzed runoff series has a high long-term memory.This is in agreement with the observed temporal behavior shown by the correlogram.In Figure 3b it can be clearly observed that the zone of the highest dependency is mainly located in the time-lag interval [0,7].
An in-depth analysis of the correlogram reveals isolated points that may provide uncertainty on time dependence (see time-lags 4, 7 and 8).In these points, the correlation coefficient (rk) is very close to the Anderson limits, rk 4 = 0.225 [−0.257, 0.226], rk 7 = 0.216 [−0.263, 0.231] and rk 8 = 0.231 versus [−0.265, 0.232].In addition, from rk 8 to rk 10 a slight trend toward an independent area of the correlogram is also observed.This is due to the average and static behavior that a correlogram offers.However, these doubts on the dependent nature of runoff series and its time-horizon will be solved by means of the dynamic analysis through the Causality.In addition, from rk8 to rk10 a slight trend toward an independent area of the correlogram is also observed.This is due to the average and static behavior that a correlogram offers.However, these doubts on the dependent nature of runoff series and its time-horizon will be solved by means of the dynamic analysis through the Causality.Synthetic series, obtained through a parsimonious and unconditioned ARMA (1, 1) model, maintain the main statistical parameters of the runoff series (Table 2).In this sense, Figure 4 shows the spectrum of generated synthetic series to populate and train the Causal Reasoning through Bayesian modelling.Synthetic series, obtained through a parsimonious and unconditioned ARMA (1,1) model, maintain the main statistical parameters of the runoff series (Table 2).In this sense, Figure 4 shows the spectrum of generated synthetic series to populate and train the Causal Reasoning through Bayesian modelling.

Analysis of the Temporal Conditionality through Causal Reasoning
The temporal dependence propagation (mitigation) for the period 1945/46 to 2013/14 was evaluated by means of a dynamic analysis of the relative percentage of change evolution over the whole time period.Figure 5 shows the discovered dependence through the maximum (MAX or

Analysis of the Temporal Conditionality through Causal Reasoning
The temporal dependence propagation (mitigation) for the period 1945/46 to 2013/14 was evaluated by means of a dynamic analysis of the relative percentage of change evolution over the whole time period.Figure 5 shows the discovered dependence through the maximum (MAX or positive) and minimum (MIN or negative) wrap-around shapes and fourth-order functions, which describe the temporal behavior, these last ones with high determination coefficients (R 2 = 0.99 in both cases).Additionally, the analytical power and suitability of Causal Reasoning applied to complex natural systems are clearly revealed in Figure 6, in which a large number of non-trivial dependency relationships (time-lag >1) amongst water years (decision variables) are revealed.Every detected connection shows an independent relationship that is measured by its p-value.Finally, from the graph of marginal dependencies and for a 100% independence level threshold and Equation 3, Figure 7 presents, in matrix form, the hydrological memory of the basin in terms of the runoff independence/dependence strength, over time.It is also worthy highlighting the lack of symmetry in the Figure 5, with a positive (MAX) wrap-around evidently dominant over the negative one (practically null).This reveals a clear dependence behavior which coincides with that observed in the correlogram (see Figure 3b).However, the dynamic propagation of dependence by Causal Reasoning encompasses a time-horizon of exclusively 7 years (from maximum value of relative percentage change (time-lag = 0 [0,1]) to minimum value (y-axis equal zero; time-lag = 7 [6,7]), versus 10 years shown by the average and static behavior that the correlogram offers (Figure 3b).Furthermore, uncertain points on time dependence shown in the correlogram (rk 4 and rk 7 ) have been eliminated by Causal Reasoning.This is shown by the dominant positive wrap-around that reveals that dependence exists but with less force.In this sense, both indeterminacies (rk 4 and rk 7 ) might be due to the same aforementioned reasons for the observed difference in the time horizon.
A detailed analysis of the graph of dependence propagation highlights a double dependence behavior.The first one, between interval [0,4], of greater dependence, presents a rapid mitigation of dependence from a relative percentage of change of 4,612.7 (time-lag = 0) to 167.7 (time-lag = 4), and where the maximum dependence is observed in the interval [0,1].In contrast, the second one, in the interval [4,7] and with a reduced relevance, shows a smooth mitigation and a non-significant increase in the final interval [6,7].These discovered, different and continuous behaviors might indicate a temporal behavior both in the short ([0,4]) and medium ( [4,7]) term that is not detected by the correlogram.
Additionally, the analytical power and suitability of Causal Reasoning applied to complex natural systems are clearly revealed in Figure 6, in which a large number of non-trivial dependency relationships (time-lag >1) amongst water years (decision variables) are revealed.Every detected connection shows an independent relationship that is measured by its p-value.
Finally, from the graph of marginal dependencies and for a 100% independence level threshold and Equation (3), Figure 7 presents, in matrix form, the hydrological memory of the basin in terms of the runoff independence/dependence strength, over time.Additionally, the analytical power and suitability of Causal Reasoning applied to complex natural systems are clearly revealed in Figure 6, in which a large number of non-trivial dependency relationships (time-lag >1) amongst water years (decision variables) are revealed.Every detected connection shows an independent relationship that is measured by its p-value.Finally, from the graph of marginal dependencies and for a 100% independence level threshold and Equation 3, Figure 7 presents, in matrix form, the hydrological memory of the basin in terms of the runoff independence/dependence strength, over time.For the color code please see Table 1.

Management Scenarios
Applying Equations 4 and 5 over runoff series, the two temporal fractions (TCR and TNCR) were determined and quantified in agreement with the graph of mitigation/propagation of temporal dependence (see Figure 5).
Figure 8 shows the different generated dynamic management scenarios.In addition, Table 3 For the color code please see Table 1.

Management Scenarios
Applying Equations ( 4) and ( 5) over runoff series, the two temporal fractions (TCR and TNCR) were determined and quantified in agreement with the graph of mitigation/propagation of temporal dependence (see Figure 5).
Figure 8 shows the different generated dynamic management scenarios.In addition, Table 3 summarizes the evolution of both temporal fractions with respect to the mean value of runoff.In each time-lag, it can be clearly seen, as the percentages of TCR exhibit high values (always above 65% and up to mean value of 86.1%), what agrees with the Hurst coefficient (0.84).Furthermore, the dual pattern of revealed behavior by Figure 5 in two intervals [0,4] and [4,7] is also observed.The discovery potential that this methodological framework offers for the management of the water resources of a basin is shown in the Tables 4 and 5. Regarding the ARMA models coefficients of the TCR fraction, their range does not show a significant variability.However, slight nuances are evident in their behavior because of the propagation of time dependence.Classical analysis does not detect this because the coefficients are unique and fixed in each model.In contrast, the TNCR fraction variability is greater, especially in the ∅ 1 coefficient (from 0.9979 to 0.9993) versus the fixed coefficients of the historical series (0.9985 and 0.9992).This would allow for a better characterization of this uninfluenced by time fraction, and therefore, the uncertainty of temporal behavior would be reduced.Figure 9 shows the design of the dynamic predictive model over the time-horizon of one (1) year (target water year 2014/15).This was performed by means of the conditioned data (Table 3), TCR and TNCR fractions time series for a one-year time dependency (Dep-1; Figure 8) and two parsimonious ARMA (1,0) models, defined each of them over those runoff temporal-fractions.It should be noted that model selection was carried out by Akaike Information Criterion, due to AIC of the ARMA models from (1.0) to (2.2) reported a slight difference in each one of the fractions (TCR: [−5.8267, −5.8265]; TNCR [−5.8469, −5.8466]; see Table 5).
approach and a detailed one developed through a definition of TNCR intervals which is an uninfluenced fraction of the time.From this approach, this adds an uncertainty component to the prediction because this fraction does not depend on time.The recorded value of target water year to assess the predictive reliability, 28.4 Hm³ was supplied by the Jucar River Basin Authority (please, see Reference [41]).On the other hand, the reliability of the predictive model is shown in Table 6, both a general approach and a detailed one developed through a definition of TNCR intervals which is an uninfluenced fraction of the time.From this approach, this adds an uncertainty component to the prediction because this fraction does not depend on time.The recorded value of target water year to assess the predictive reliability, 28.4 Hm 3 was supplied by the Jucar River Basin Authority (please, see Reference [41]).

Discussion
This methodology allows knowing, in detail and jointly, the behavior of water resources in the short and medium term (Figure 5; [0.4] and [4.7]), as well as their time-horizon (7 years versus 10 of the correlogram; Figures 3b and 5).In this sense, it is worth highlighting the coherence between the Hurst coefficient (0.84) and the results obtained (asymmetry in Figure 5).It is also remarkable the existing agreement between the dual pattern of behavior presented in Figure 5 (intervals [0,4] and [4,7]) and the analysis of the TCR and TNCR fractions´temporal behavior evolution shown in Table 3 (time-lags 1 to 4; 84.8% to 70.1%).
Furthermore, the potential of this methodological framework for the basin water resources management largely became clear when the temporal behavior nuances of each fraction series (Tables 4  and 5) were discovered.It is important to mention that these nuances cannot be detected by the classical approaches.This was revealed by the variability of the ARMA model coefficients, especially noteworthy in the case of TNCR Fraction (Table 5).
Regarding the predictive model, it shows a high level of reliability, 82.4% in general perspective, and within 85% to 90% through a detailed approach ([20.11, 30.23] and [22.26, 34.50] for a time-dependency of one-year; see Table 6).
Moreover, detailed analysis over the temporal fractions presents crucial and convenient implications over dimensioning, and control operations of reservoirs and dams.This is clearly shown on the uncertainty that TNCR fraction adds to the predictive model.In a sense, for a runoff of 28.4 Hm 3 in the target water year, there is an uncertainty among 5.06 and 6.12 Hm 2 due to the fraction not influenced by time, approximately 17.8 or 21.5% of the total runoff, respectively.Therefore, in a management context, that volume should be jointly provided by the reservoir-dam and its control operations in the previous year.In this way, the failure capacity of the infrastructure is reduced or minimized.
Although reservoir dimensioning involves knowing multiple key aspects such as capacity-areaelevation curves, downstream impacts and sediments inflow, among others, this methodological framework might effectively contribute to achieving a better dimensioning of reservoirs and dams.This will necessarily improve service guarantees through detailed knowledge of the runoff temporal-fractions (TCR and TNCR), especially TNCR fraction, which it is not temporally dependent.
Furthermore, this methodology opens new perspectives for building dynamic and stochastic hydrological predictive models.Also, more river basins sustainable management approaches are expected to be designed.This may be produced through more reliable control operation of reservoirs within the current context of global warming, and for high temporal dependent river basins.
Future research will incorporate the analysis of spatial dimension behavior in the design of causal models.

Conclusions
This research supposes a breakthrough in the temporal characterization of runoff series.Here, Causality has discovered a hidden and logical structure of temporal interdependence that inherently underlie the hydrological series.This represents a latent behavior pattern that was waiting to be discovered and that the classical approaches had not been able to reveal.This was performed by means of the coupling of a stochastic hydrological parametric, parsimonious and unconditioned ARMA model, and Causal Reasoning.The latest was supported by Bayesian modelling, which is a powerful AI technique based on probabilistic graphical model.
The methodological framework makes possible to generate dynamic management scenarios according to the propagation (mitigation/evolution) of the dependence (time-lags), taking into account the runoff dependence strength over time.In that sense, two opposite temporal-fractions, the Temporally Conditioned Runoff (TCR) fraction, and Temporally Non-Conditioned Runoff (TNCR) fraction were determined and quantified within the runoff series.Furthermore, the observed behavior trends of TCR/TNCR fractions made it possible to build a predictive runoff model with a high level of reliability.
The application of this new analytical strategy over annual runoff series applied to dimensioning of reservoirs, dams, and control operation of reservoirs seems quite straightforward.In this sense, a design more adjusted to the real needs of reservoirs-dams through a better knowledge and prediction on temporal behavior of TCR and TNCR fractions is expected to be achieved.This will lead to a reconsideration of the capacity of the storage reservoirs, especially their convenience in the current context of global warming and for high temporal dependent river basins.This would imply substantial economic and environmental savings in the future.
This research has demonstrated that past information provides prior knowledge of the future with a high degree of reliability.Furthermore, this research reinforces the suitability of Causality (Causal Reasoning) in modelling the temporal behavior of the water resources of a highly dependent basin from a dynamic and stochastic perspective.

Figure 1 .
Figure 1.(a) Situation map; (b) Jucar River Basin map and case study situation; (c) Historical annual series of runoff (time period: 1945/46 to 2013/14).Note: Time series are referred according to water year in Spain where the first month is October and the last is September of following year.

Figure 3 .
Figure 3. (a) Main statistical parameters of runoff series.(b) Correlogram and Anderson probability limits for an independent series (95 percent probability level).

Figure 3 .
Figure 3. (a) Main statistical parameters of runoff series.(b) Correlogram and Anderson probability limits for an independent series (95 percent probability level).

Figure 4 .
Figure 4. Spectrum of synthetic series to populate the Causal Reasoning versus historical runoff.

Figure 4 .
Figure 4. Spectrum of synthetic series to populate the Causal Reasoning versus historical runoff.

Water 2019 , 19 Figure 5 .
Figure 5. Causal modelling of the mitigation/propagation of temporal dependence by an analysis of the evolution of the relative percentage of change over the time period 1945/46 to 2013/14.

Figure 6 .
Figure 6.Graph of Marginal Dependencies.Note: The displayed threshold of level of independence is 0.05 (up to 95% of dependency relationships are shown amongst water years or decision variables).

Figure 5 .
Figure 5. Causal modelling of the mitigation/propagation of temporal dependence by an analysis of the evolution of the relative percentage of change over the time period 1945/46 to 2013/14.
the evolution of the relative percentage of change over the time period 1945/46 to 2013/14.

Figure 6 .
Figure 6.Graph of Marginal Dependencies.Note: The displayed threshold of level of independence is 0.05 (up to 95% of dependency relationships are shown amongst water years or decision variables).

Figure 6 .
Figure 6.Graph of Marginal Dependencies.Note: The displayed threshold of level of independence is 0.05 (up to 95% of dependency relationships are shown amongst water years or decision variables).Water 2019, 11, x FOR PEER REVIEW 11 of 19

Figure 7 .
Figure 7. (a) Matrix of independence.The cell value is p-value.(b) Matrix of dependence.The cells show dependence rate (DR) value of each relationship.Row values indicate year interactions with respect to target year (main diagonal cell).Note: In order to achieve a better knowledge about the results obtained a part of the matrixes is only shown (1945/46 to 1958/59), nevertheless the results encompass the whole time period of the runoff series.For the color code please see Table 1.

Figure 7 .
Figure 7. (a) Matrix of independence.The cell value is p-value.(b) Matrix of dependence.The cells show dependence rate (DR) value of each relationship.Row values indicate year interactions with respect to target year (main diagonal cell).Note: In order to achieve a better knowledge about the results obtained a part of the matrixes is only shown (1945/46 to 1958/59), nevertheless the results encompass the whole time period of the runoff series.For the color code please see Table 1.

Figure 8 .
Figure 8. Dynamic management scenarios.Please note Dependence Propagation (3 year) is referred time-lag equal to 3.

Figure 8 .
Figure 8. Dynamic management scenarios.Please note Dependence Propagation (3 year) is referred time-lag equal to 3.

Figure 9 .
Figure 9. Preliminary implementation of the predictive model for a Dependence Propagation one year (Dep-1).

Figure 9 .
Figure 9. Preliminary implementation of the predictive model for a Dependence Propagation one year (Dep-1).

Table 1 .
Dependence Rate (DR).Definition of classes and intervals.

Historical runoff series Average of all annual synthetic series
1 Classic skewness coefficient.

Historical Runoff Series Average of All Annual Synthetic Series
1 Classic skewness coefficient.Water 2019, 11, x FOR PEER REVIEW 9 of 19

Table 3 .
Temporally Conditioned Runoff (TCR) and Temporally Non-conditioned Runoff (TNCR) Fractions.Analysis of temporal behavior for a reference threshold 2 .
2Reference threshold: Mean value of the runoff.

Table 4 .
TCR Fraction.Evolution of the main statistical parameters according to time dependence.

Table 5 .
ARMA model selection.parameters and the coefficients.

Table 6 .
Prediction Analysis for a Dependence Propagation one year (Dep-1).