New Insights on Land Surface-Atmosphere Feedbacks over Tropical South America at Interannual Timescales

We present a simplified overview of land-atmosphere feedbacks at interannual timescales over tropical South America as structural sets of linkages among surface air temperature (T), specific humidity at 925 hPa (q925), volumetric soil water content (Θ), precipitation (P), and evaporation (E), at monthly scale during 1979–2010. Applying a Maximum Covariance Analysis (MCA), we identify the modes of greatest interannual covariability in the datasets. Time series extracted from the MCAs were used to quantify linear and non-linear metrics at up to six-month lags to establish connections among variables. All sets of metrics were summarized as graphs (Graph Theory) grouped according to their highest ENSO-degree association. The core of ENSO-activated interactions is located in the Amazon River basin and in the Magdalena-Cauca River basin in Colombia. Within the identified multivariate structure, Θ enhances the interannual connectivity since it often exhibits two-way feedbacks with the whole set of variables. That is, Θ is a key variable in defining the spatiotemporal patterns of P and E at interannual time-scales. For both the simultaneous and lagged analysis, T activates non-linear associations with q925 and Θ. Under the ENSO influence, T is a key variable to diagnose the dynamics of interannual feedbacks of the lower troposphere and soil interfaces over tropical South America. ENSO increases the interannual connectivity and memory of the feedback mechanisms.


Introduction
The humid tropics receive a large amount of net radiation and water vapor, and consequently, intense heat and humidity fluxes dominate the interactions between soil and the lower atmosphere [1]. The excess of net radiation is balanced through latent and sensible heat fluxes [2][3][4]. By controlling the partition of these fluxes, soil moisture modulates diverse land-atmosphere feedbacks (LAFs); similarly, the rates of change between dry and wet conditions in the soil necessarily impact surface temperatures [5,6].
The study of LAFs has been approached from physical models and numerical experiments [7][8][9][10][11][12][13]; analysis of observations and models with statistical tools [3,[14][15][16][17]; and traces of moisture trajectories [18,19], among others. An important body of literature has focused on the role of vegetation and land uses in the dynamics of LAFs [20][21][22][23][24], and the conditions under which LAFs determine the stability of the lower atmosphere [25][26][27]. According to Figure 1, soil moisture regulates the entire LAF system through its direct linkage to evapotranspiration and precipitation [30,49,50]. Besides, soil moisture also influences soil temperature through albedo and latent heat flux [51], and at the same time non-linear vegetation-induced turbulence, determined by vegetation stature and cover fraction, significantly influences ET-soil moisture relationships [52]. Another important feedback is the influence of soil temperature, in connection with air temperature, sensible heat flux, and evaporation [48].
The scheme of relationships between state and process variables shown in Figure 1 is a hierarchical adaptation using Graph Theory of the original figure presented by Brubaker and Entekhabi [48]). Graph Theory is a specific branch of mathematics that studies the interconnections between objects. An interpretation of couplings, feedbacks, and interactions via Graph Theory to further understand LAFs was proposed in previous studies [42]. According to such interpretation, the concept of feedback with graphs coincides with the type of bidirectional link established in Figure 1 between evaporation and soil moisture; on the contrary, the precipitation-soil moisture link is defined as a coupling or a hierarchical interaction whereby one variable dominates over the other. The structure of links shown in Figure 1, derived from graph theory, provides a clearer understanding of the processes and mechanisms involved in LAFs [53][54][55][56].
The aim of this paper is to explore new ideas and tools to advance our understanding of land surface-atmospheric feedbacks in TropSA at interannual timescales. To achieve this goal, we propose to integrate classical analyses of climatic fields with more recent methods of information theory and Graph Theory. Our approach involves linear and non-linear analysis, including metrics such as causality between variables. An explicit goal of this paper is to study the essential role of such feedbacks to understand the mechanisms involved in water and heat anomalies in TropSA associated with El Niño-Southern Oscillation (ENSO) [44,57].
The study is organized as follows. Data and methods are presented in Section 2. In Section 3, we define spatiotemporal patterns of maximum covariance between pairwise selected variables over TropSA during the period 1979-2010, and, in order to estimate the interannual linear and non-linear connectivity of the set of variables, we use metrics such as correlation and causalities among the representative time series of the estimated maximum covariance patterns (Sections 3.2 and 3.3). Using the same time series, we select the patterns with the highest influence of ENSO. Finally, with these patterns, we evaluate the structure of relations variables using elements of graph theory (Section 3.3).

Data
We focus our study on tropical South America (15 • N-20 • S and 82 • W-40 • W, Figure 2) taking some of the variables proposed by Brubaker and Entekhabi [48] to the study of LAFs. Our analysis involves the interactions among surface air temperature (T), specific humidity at 925 hPa (q 925 ), and soil moisture (Θ) (as state variables), precipitation (P) and evaporation (E) (as process variables) [7].
Precipitation comes from the Global Precipitation (GPCC) Full Data Product version 7 (ftp: //ftp.dwd.de/pub/data/gpcc/html/fulldata_v7_doi_download.html) [58]; the GPCC monthly precipitation fields are suitable to study hydrological and atmospheric linkages at regional scales on land. From the ERA-Interim Reanalysis [59], we use monthly fields of T, q 925 , and E, the latter represented by the Instantaneous Moisture Flux. Volumetric Soil Water content (Θ) is obtained from the ERA-Interim/Land Reanalysis V.2 [60,61]. This latter product is focused on those land processes not well represented in the ERA-Interim Reanalysis and is just available for the period 1979-2010 [62], which justifies this time span for our study. The monthly climatic fields were originally downloaded at 1 • × 1 • spatial resolution.
To extract the interannual variability, we removed the annual cycle from all data sets by a simple standardization procedure (subtracting the monthly mean and scaling by the monthly standard deviation). Data from the ERA-Interim reanalysis require a split climatology associated with systematic discontinuities in the water balance due to changes in the sensor channel of the Microwave Special Sensor (SSM/I) [63]. Thus, we estimated separate means and standard deviations for the periods 1979-1991, 1992-2002, and 2003-2010 to estimate monthly anomalies.
We also use the Niño 3.4 index to represent the dynamics of ENSO over the tropical Pacific Ocean at interannual timescales, and monthly fields of SST from the Hadley Center (HADISST) to map the statistically significant correlations between hydro-climatic fields over TropSA.

Maximum Covariance Analysis
Since we aim at studying the spatiotemporal covariance between all possible pairs of variables, 10 possible combinations are possible. The covariance represents the interaction between each pair of variables at interannual timescales. We use the matrices of standardized data X and Y during the 384 common months in the period 1979-2010. Then, we estimate the covariance matrix, C XY , as: C XY was factorized employing a singular value decomposition (SVD) to extract those patterns explaining the maximum covariance fraction between X and Y at the interannual timescale. This method is also known as Maximum Covariance Analysis (MCA) [54,64,65]. Through the MCA, C XY can be decomposed as: In this case, the first singular factor, u 1 σ 1 v T 1 , is the leading MCA-mode of covariance between X and Y. The magnitude of the first singular value, σ 1 , indicates the amount of covariance between both variables, which is captured by the leading mode. The second singular factor, u 2 σ 2 v T 2 , is orthogonal to the first one, and maximizes the remaining portion of the covariance, and so on until explaining all the covariance contained in C XY . A measure of the degree to which both variables maximize their temporal covariance within each mode, is provided through the estimation of the MCA-series x k and y k , by projecting the matrices u k and v k (containing the singular vectors of C XY ) on the original matrices X and Y (Equation (3)), as: Each mode contains a fraction of the square covariance, f XY k , defined as: f XY k represents the magnitude of the covariance that is explained by each mode between each pair of variables. The sum of the k relative fractions is denoted by F XY k . Finally, we estimate the spatial patterns associated with X and Y for each MCA-mode. To that end, we use correlation maps between the series x k and y k and the matrices of standardized data X and Y.

A Graph Model and Feedbacks at Interannual Timescales
The combined usage of MCA and Graph Theory has been recently implemented in some studies [53,54,56]. Graph models have also been used to evaluate climatic interactions in a context of causality [53,56]. For our purposes, we selected the interaction of the MCA-series according to the interannual mode that determines their maximum covariance. The criterion to group the MCA-series of each matrix is defined by the highest strength of correlation with the Niño 3.4 index.
For lags from 1 to 6 months, a graph model is constructed to study the connectivity at the interannual timescale for each MCA-mode between variables. We use the concept of links between two variables x and y [49] to construct the graph model, Γ, as: where the nodes V are the variables with maximum covariance under the highest ENSO influence. The weight of the edges, E, are estimates of the correlation and causality metrics between the x and y from each MCA-mode (see Supplementary Materials). The graph depicting the linear associations is constructed by assigning the edge weights as correlations (ρ). A graph of non-linear associations is constructed with causalities (τ) (see Supplementary Materials).
Additionally, we introduce a metric to quantify the ratio of bi-directional causalities between two series y and x (Equation (6)), as: A value of Φ y←→x = 1 indicates a bidirectional feedback, while Φ y←→x = 1 quantifies the degree of coupling between both variables. For our purposes, we established a threshold τ y←→x > |1%| to define the relative causalities quantifying the interaction strength. With this restriction, we use only the strongest causalities in each pair of relationships.
Estimation of simultaneous correlation does not provide information about directionality, and thus the resulting linking graph is called an Undirected Graph. Otherwise, by applying bidirectional metrics, such as causality described in Supplementary Materials [66], it is possible to establish a Directed Graph, denoting directions with arrows ( Figure 3). Lagged correlations and causalities also quantify bidirectional characteristics. Figure 3 shows a basic relationship between two-time series x k and y k as a schematic view in two-time steps: x(t − 1) influences y(t) and y(t − 1) influences x(t). The red arrow in Figure 3 connects the series when x leads y, and the blue arrow when y leads x. Hence, in both cases, we systematically estimate correlations, ρ(x(t), y(t − 1)) and ρ(x(t − 1), y(t)), and causalities, τ (x(t),y(t−1)) and τ (x(t−1),y(t)) .
We define an edge magnitude threshold to construct correlation and causality graphs to focus on the strongest connections among the studied variables on the highest ENSO-related covariance modes. In all cases, we selected the percentile of 50% of the empirical probability distribution of causalities and correlations values for each time lag.
The adjacency matrix was extracted to study the links structure among variables in each graph. Factorizing each adjacency matrix by applying an SVD and following the HITS algorithm proposed by [67], we systematically extracted the highest entries of the leading singular vectors pair from each matrix factorization. In this context, the adjacency matrix directly connects the whole set of graph nodes crosslinking their columns (emitters) and rows (receivers). As the adjacency matrix defines all the possible associations among nodes, the variable (node) corresponding to the maximum vector u 1 value (emitter) transfers the highest amount of information to the whole set of variables, while the maximum vector v 1 value is the best information receiver.

Maximum Covariance Analysis among Variables and Strength of ENSO Influence
The three main MCA-modes are able to capture a large portion of the total covariance among the interactions of the five variables in the study region (at least 76%). Table 1 shows results hierarchically sorted according to the covariance fraction from the leading MCA-mode of each pair. The T-q 925 interaction presents the highest fraction in the leading mode corresponding to 73% of its total covariance. This is because both variables are strongly linked by the Clausius-Clapeyron formulation between temperature and saturation water vapor pressure, especially in the humid tropics [1]. Table 1. Three leading MCA-modes of all possible covariance matrices C XY . Fraction of relative ( f XY k ), cumulative square covariance (F XY k ), and simultaneous correlation between MCA-series, ρ(x k , y k ). Pairs of variables are sorted according with the value of the first fraction. The mechanisms of the hydrological cycle (P-Θ-E) are concentrated in the main two interannual modes over the study region collecting up to 71% of the covariance. Compared with the whole set of interactions, E-P and P-Θ concentrate the lowest fraction in the leading MCA-mode and the highest fraction in the second one (Table 1). However, the interaction of E and Θ share more covariance than P-Θ in the leading mode. Such a strong covariability between soil moisture (Θ) and evaporation (E) reinforces the idea that droughts are also associated with LAFs at local scales, and with teleconnections at larger scales in TropSA [23,68,69].
We also assessed the influence of ENSO on the covariance modes correlating the MCA-series with the Niño 3.4 index ( Table 2). According to the magnitude of these correlations, we identify the largest ENSO-influenced mode at each pair of variables (in bold). The interaction Θ-T shows the highest correlation with N3.4 among the paired variables (ρ = 0.58) in the leading MCA-mode. Soil moisture (Θ) influences the atmospheric response to surface warming (T), and, in combination with atmospheric moisture (q 925 ), Θ contributes to determining the spatiotemporal pattern of P [70]. Regularly, those pairs involving T exhibit the highest correlation values with the Niño 3.4 index ( Table 2); this means that surface temperature is essential to propagate the ENSO-effect on studied variables. Table 2. Correlation between MCA-series with the Niño 3.4 index, ρ(x k , N3.4) and ρ(y k , N3.4) for the three leading MCA-modes. Highest correlation value of each pair is denoted in bold.
Each C XY including P (P-q 925 , P-T, E-P, and P-Θ) consistently shows the highest correlations with N3.4 in the second MCA-mode ( Table 2). The P-q 925 pair shows a clear difference of correlation with N3.4 among the three modes; in that case, the second mode of P-q 925 is best correlated with N3. 4. Figure 4b,c shows the typical spatial pattern of the highest-related ENSO modes corresponding to the interaction between evaporation (E) and soil moisture (Θ). The core of this pattern is located in the Amazon River basin and is connected to the Magdalena-Cauca River basin and the western Pacific coast of Colombia where the Choco Jet has a strong effect on the hydroclimate variability [68] ( Figure 4B,C). Some regions of the Orinoco and Tocantins River basins do not show statistically significant values (p < 0.01) under this mode. Spatially, the pattern of correlations suggests that the influences of E on Θ are higher than in the reverse case (Figure 4b,c). Both maps exhibit negative correlations, suggesting that negative (positive) anomalies of Θ correspond to negative (positive) anomalies in E. Temporally, values of one-month lag correlation coefficients suggest that the effect of Θ on E is higher than in the reverse sense (Figure 4d). This effect is coherent but less evident at two-month lags. The cross-correlation between both series is considerable up to six-month lags.
Complementary, Figure 5 shows the global SST-correlation map that corresponds to the high ENSO-related using the MCA-series presented in Figure 4a. The SST-pattern exhibits highly statistically significant correlations with the Niño regions acting in phase (same sign of coefficients) with the Tropical North Atlantic (TNA), and with the Indian Ocean SSTs (IO) ( Figure 5). The pattern also exhibits the influence of the Pacific Decadal Oscillation (PDO) acting in phase-shift with Niño regions correlations (reverse sign of coefficients). This conjoint effect is triggered because the Amazon River basin occupies the largest portion of TropSA and several SST-modes determine the characteristics and magnitude of the interannual anomalies therein [57,[71][72][73].

Non-Linear Analysis of the MCA-Series and Feedbacks
We estimate the relative causalities τ y→x using all pairs of MCA-series with the highest influence of ENSO. Table 3 presents the non-linear metrics among those MCA-series. The interaction q 925 -T shows the largest relative causalities from the set of estimated values. Specifically, T reaches up to −23.3% of causality on q 925 . The negative value indicates that T favors the predictability of q 925 under the influence of ENSO [70]. After q 925 -T, the P-Θ interaction shows the largest causality value. Θ consistently subordinates P (τ y→x = 21.6%), while the reverse causality is 15.2%, implying a coupling Φ VSW←→PRC = 1.43. These non-linear connections suggest that Θ is fundamental in modulating the interannual variability of P over the studied region under the ENSO influence. In Table 3, we also highlight in bold those values of Φ (Equation (6)) within the interval [0.9, 1.1] defined as a bidirectional feedback. These values are more frequent under the highest ENSO-influenced mode. Surface temperature (T) acts in two of those feedbacks (q 925 -T, and T-Θ). Additionally, the interaction T-Θ shows a perfect bidirectional feedback (Table 3). T is essential in capturing the non-linear interactions at interannual timescale within the studied variables in the simultaneous analysis.

Visualizing Linear and Non-Linear Dependences Using Graph Theory
For our study, Graph Theory has the advantage of representing connections among the evaluated variables, assuming that they represent the LAF feedbacks in Tropical South America. In this sense, such feedbacks can be understood as the degree to which the atmosphere responds to anomalies in the land surface state [15]. Soil moisture has been identified as a key variable in several studies and especially in the terrestrial segment of the hydrological cycle i.e., the evaporation-soil moisture interaction [4,[8][9][10]15,16]. Using Graph Theory, we verify the interdependences of several variables in the context of these complex interactions at interannual timescales. Table 4 presents the median of correlation and relative causality at each lag for all possible interactions among the five selected variables. Based on these thresholds, we construct graphs of linear ( Figure 6) and non-linear (Figure 7) linkages. To construct the graphs, we only used the values higher than those thresholds. The couple of red nodes in each graph indicate the variables with the highest value in the leading singular vectors u 1 and v 1 estimated from the SVD applied to the graph adjacency matrix. We also show those variables as red horizontal bars in Figures 6 and 7. Red paired-edges in each graph represent correlations and causality ratios (Equation (6)) ranging between 0.9 and 1.1. Table 4. Threshold of edges magnitude to construct correlation and causality graphs (percentile of 50% of the empirical probability distribution of causalities and correlations per lag).  Figure 6 shows the linear graph constructed with linkages at one-, two-and three-month lags under the highest ENSO-related MCA modes. The interactions among the five variables are connected with the spatial patterns showed in Figure 4b,c over TropSA and the SST-pattern established in Figure 5. Linear interactions among the variables depict that soil moisture (Θ) structurally regulates the mechanisms of interdependences among the five variables for antecedent/subsequent times [9].

Metric Type
Specifically, Figure 6a shows the graphs whose edges are one-month lag correlations higher than the median (ρ = 0.43, Table 4). The structural interaction between soil moisture-atmospheric humidity (Θ and q 925 ) corresponds to the highest entries of leading singular vectors -u 1 and v 1 -of the graph matrix (Figure 6a). The high degree of response of the atmosphere to anomalies in the land surface state is taking into account all the interdependences among the five selected variables. Our results point out that ENSO favors an important chain of linear feedbacks among T, q 925 , and E. This means that land-atmosphere interactions are fundamental driving mechanisms of interannual variability and delayed effects on hydrologic response over the study region [42,45]. Led by the ENSO-forcing, Θ presents the highest correlation value with P among the variables (ρ = 0.69, Figure 6a) implying that the interannual anomalies of soil moisture are strongly related with the subsequent anomalies at one month of precipitation over the region [42]. Surface temperature anomalies are also directly regulated by soil moisture anomalies and through linear feedbacks with E and q 925 (Figure 6a). From this same graph, we infer that antecedent anomalies of q 925 do not influence directly the subsequent anomalies of Θ but through changes in T; instead, one-month lagged anomalies of Θ linearly controls the anomalies in q 925 .  Figure 6b shows the graph constructed with two-month lags correlations also associated with the highest ENSO-influenced mode. Again, this graph reflects a relational structure dominated by Θ as the best emitter (maximum value of u 1 ); but, in this case, evaporation (E) is the best receiver (maximum value of v 1 ). For two-month lags, T establishes linear feedbacks to connect q 925 with Θ. Again, the highest correlation of the entire set shows the predominance of Θ on P (ρ = 0.53). Soil moisture anomalies regulate the interdependence of P anomalies with T, q 925, and E anomalies at interannual timescale.
Finally, Figure 6c shows the graph corresponding to the three-month lags revealing a structural relation between surface temperature and humidity at 925 hPa (T-q 925 ). As shown by one-and two-month lags, the structure of the graph represented in a hierarchical form demonstrates that Θ constitutes a heterogeneous node connecting the anomalies of E, q 925, and T with the subsequent anomalies of precipitation (Figure 6c). However, in this case the correlations P-Θ are not the highest among the variables, and thus the structure of the graph is inverted with respect to Figure 6a,b. The linear feedbacks are identified for the pairs E-Θ, Θ-q 925 and q 925 -T (Figure 6c). Among the studied variables, T and Θ exhibit the highest amount of linear feedbacks, including the correlation among them at one and two month-lag, and with E. Figure 7 presents graphs whose edges denote causality relations (τ y←→x , non-linear associations). Figure 7a shows the simultaneous (non-lagged) graph for the causalities reported in Table 3. Contrary to correlations, relative causalities allow us to examine the relations structure without lagging. Because of non-linearity, relative causalities depend on direction between the two variables, given by τ y→x and τ x→y (Table 4). Surprisingly, the linkages structure is similar to that depicted by linear graphs of one-month lag (Figure 6a) and two month-lags (Figure 6b). In this case, Θ also works as a transferring node variable connecting P with E, q 925, and T (Figure 7a). The highest causality of this graph of concurrent interactions corresponds to the T → q 925 relationship (23.3%). In addition, the structure of this graph is nearly equally dominated by Θ as the best emitter (u 1 ), while the best receptor is q 925 (Figure 7a). This means that Θ dominates the non-linear relationships with P, q 925, and T at interannual timescales and without lagging. It is remarkable that under this structure of causalities, evaporation (E) is the only variable which non-linearly controls Θ. This agrees with the essential role of evaporation in the surface energy partitioning in wetter environments [74], taking into account that under this non-linear framework, precipitation presents a relative causality on soil moisture (Figure 7a). In general, Θ is essential over Tropical South America to structurally feedback on q 925 under the highest influence of ENSO at concurrent series (non-linear links).

Non-Linear Graphs
While the simultaneous causalities are evidence that Θ dominates q 925 (Figure 7a), T exerts a non-linear control on Θ at one-and two-month lags (Figure 7b,c). In this way, non-linear associations reveal that T is essentially strengthening the identified feedbacks under the ENSO effect [6]. For both the simultaneous and lagged analyses, surface temperature (T) exhibits non-linear associations with the atmospheric moisture (q 925 ) and soil moisture (Θ) causality graphs. Under the ENSO influence, T is not only a key variable in the diagnostics of the dynamics of interannual LAFs but also has a substantial role in the dynamics and thermodynamics of the lower troposphere and soil interfaces over TropSA [2,3,40,44,69,75]. Besides, interannual anomalies in P are associated with perturbations of the surface water balance at the regional scale by controlling T interactions with Θ and E [10,17,30,31,45,49,[76][77][78].
Comparing the graphs under one-month lag for linear ( Figure 6a) and non-linear (Figure 7b) connections, we found different relational structures. Under the non-linear structure, Θ receives relative causalities from the rest of variables (P, E, T, and q 925 ) with a highest value coming from temperature (τ T→Θ = 11.4%, Figure 7b). Another important fact is the causality that precipitation has on temperature (τ P→T = 5.45%). This implies an important non-linear connection of the antecedent anomalies of temperature and precipitation on the one-month subsequent anomalies of soil moisture.
This also contributes to understanding that the linear (1-month lag, Figure 6a) and non-linear (simultaneous, Figure 7a) structure of relationships are connected over TropSA.

Discussion
The structure of relationships illustrated by the graphs in Figure 6a-c shows the central role of soil moisture as a transitional variable between precipitation and surface processes at interannual time-scales [4,9,15]. Evaporation anomalies depend on surface conditions and these are relatively homogeneous in the region as shown by the spatial pattern depicted in Figure 4b,c [69,74]. This is the most representative ENSO-influenced area over tropical South America for all the variables. The physical mechanisms of the LAFs over Tropical South America at interannual time scales for oneto three lag-months are clearly modulated by the dynamics of soil moisture to produce the anomalies of precipitation. In this way, we reveal that the interannual anomalies witnessed in precipitation are dominated by the antecedent soil moisture anomalies for one and two-month lags. Actually, soil moisture exerts the highest influence on the studied variables and even more so on precipitation. The physical mechanisms are related to the control that soil moisture exerts over the subsequent anomalies of evaporation, atmospheric humidity and surface temperature.
One notable aspect is the fact that, under the defined thresholds presented in Table 4 for correlations, none of the variables influence rainfall as soil moisture does ( Figure 6). Only the one-month antecedent anomalies of temperature have influences on the subsequent soil moisture anomalies (Figure 6a). The role of temperature in this case is related to the influence of ENSO. The one-month lag anomalies of soil moisture dominates the subsequent anomalies of precipitation, evaporation, temperature and atmospheric humidity. This fact confirms the key role of soil moisture anomalies under the ENSO influence to define the state and dynamics of those variables (P, E, T, q 925 ) in Tropical South America at an interannual time scale.
Land surface-atmosphere feedbacks (LAFs) have been previously studied in the Amazon River basin, although their connection with the rest of tropical South America has been largely overlooked [13,20,23,25,28,34]. Here, we advance new ideas about such linkages. Besides, these LAFs have been studied for annual and intra-annual timescales, but the interannual timescale has also been largely overlooked. Most of the research at the interannual timescale has been focused on the study of large-scale atmospheric teleconnections or anomalies in the transport of moisture by winds [2,31,33,62], but the role of LAFs has not been thoroughly explored. We explored land surface-atmosphere feedbacks at interannual timescales using mainly a simplified overview of connectivity of the most relevant variables that describe these interactions as well as identifying the regions of TropSA where they are best characterized. We demonstrated that within the identified multivariate structure, Θ enhances the interannual connectivity since it exhibits two-way interactions with the whole set of variables at simultaneous and lagged times.
Results show concurrent and lagged interdependencies between soil-atmospheric humidity and precipitation over TropSA. The simultaneous interdependencies identified using relative causalities under the ENSO influence indicate that they are modulated through surface temperature (Figure 7a). This result is very interesting because it suggests a dynamic that is not possible to expose using linear analysis (correlation). The structural interaction between soil moisture (Θ) and atmospheric humidity (q 925 ) is mediated through non-linear interdependences of the surface temperature (Θ-T and T-q 925 ) as defined by Equation (6) (bi-directional causalities). At the same time, the anomalies of precipitation are highly dependent on soil moisture anomalies (Figure 7a). Results from the linear analysis show the fundamental role of soil moisture as the best emitter within the multivariate dynamics of LAFs over TropSA ( Figure 6). This role is not completely demonstrable in the non-linear graphs where temperature takes more relevance (Figure 7). The greater differences between linear and non-linear analysis correspond to the one-month lag graphs (Figures 6a and 7b). However, results indicate a connection between the structure of the linear one-month lagged and non-linear simultaneous analyses (Figures 6a and 7a). In summary, we show that linear and non-linear structural relationships are complementary to capture the dynamics of LAFs over tropical South America. We defined and estimated metrics of linear and non-linear feedbacks among several variables in a common interannual mode mainly dominated by ENSO.
Our approach allows us to quantify the role of antecedent anomalies on the related anomalies using the two-way interactions identified within the LAF framework. We present a novel approach: defining a structural scheme of relationships among the variables (via Graph Theory), using the spectral modes on the adjacency matrices of the graphs at different time lags. In addition, we introduce new metrics to quantify the ratio of bi-directional causalities between two variables (Equation (6)).
The use of MCA is useful at interannual time scales because it allows us to identify the main co-variability modes evaluating all possible relationships among the selected variables. Based on cross-correlation functions and causality, Graph Theory is a good complement in the context of land atmosphere interaction mechanisms between different processes, mainly with the purpose of visualizing and understanding the bidirectional relationships among the selected variables and also to compare linear and non-linear perspectives [53,54]. The structural relationships depicted through the use of Graph Theory are representative of the spatial pattern related to each co-variability mode.
MCA allows for simultaneous identification of main-covariability modes while also evaluating all relationships among the selected variables. This approach is important because it enables exploring of the multiple interaction and interdependencies between variables. In the context of TropSA, the MCA allows the identification of the regions and the main macroclimatic patterns influencing the interactions among the variables through the land atmosphere feedbacks. The use of network analysis to complement eigen-techniques has proved recently to be very useful in the study of climatological processes [54]. Bracco et al. [79] also discussed the use of data mining as a complement to traditional tools for climatological analysis and presented an application of Graph Theory to highlight the interconnectedness of global SST patterns among different domains. In this sense, we tried to demonstrate that, based on cross-correlation functions and causality, Graph Theory is a very good fit to complement the study of land atmosphere interactions mechanisms among different processes and state variables. In our study, we mainly used it with the purpose of visualizing and understanding the bi-directional relationships among the selected variables and also to compare linear and non-linear perspectives. The structural relationships depicted through the use of Graph Theory are representative of the spatial pattern related to each co-variability mode.
We noticed that the use of graphs is a growing and new field of Climatology [54,79,80]. In spite of the demonstrated capabilities of the approach used in the present study to reveal the temporal dynamics of two-way interactions between fundamental hydrometeorological variables, some disadvantages involve the data requirements, the capacity to estimate uncertainties, and the fact that results can be hard to interpret in a mechanistic way.

Conclusions
A multivariate analysis allowed us to further understand the role of land surface-atmosphere feedbacks over Tropical South America at interannual timescales during the period 1979 to 2010. We used elements from Linear Algebra, Information Theory and Graph Theory to infer structural relations among the most relevant processes involved in those feedbacks. In so doing, we defined the spatiotemporal patterns by a pair-wise categorization of variables through Maximum Covariance Analysis, with the aim of reducing the dimensionality of the problem, as well as to identify and quantify the most salient factors associated with ENSO. With such patterns, we evaluated the relational structure between variables using Graph Theory to estimate the linear and nonlinear connectivity at interannual timescales.
Using the ratio of linear (correlations) and nonlinear (causalities) metrics between two variables (Φ, Equation (6)), we defined a bi-directional relationship (feedback) whenever 0.9 ≤ Φ ≤ 1.1. Using such a definition, based on concepts defined in previous studies [49], we found that the causality graphs corresponding to the highest association with ENSO show several feedbacks among the studied variables for all lags between one and six months. This effect is triggered by the highest association with ENSO. Correlation metrics reveal a consistent bi-directional relation between Θ-T at one and two month-lags. Use of non-linear metrics detects the feedback Θ-T at non-lagged analysis.
The identified LAFs are enhanced (in-phase) by SSTs over the tropical Pacific and Atlantic, but also over the Indian Ocean. The regulatory action of ENSO creates a stronger memory and resilience of interactions among variables. The spectral analysis of the composite graphs, identified to study the linkages between variables, allowed us to conclude that soil moisture is a key variable in defining the spatiotemporal patterns of P and E in TropSA at interannual timescales, and as such plays a major role in regulating LAFs at interannual timescales, mainly controlled by ENSO. Although deforestation is a relevant aspect of land surface-atmosphere interactions, this topic goes much beyond the scope of our study.
Finally, within the context of the humid tropics, the identified LAFs play a major role in controlling the connection between the soil and atmosphere subsystems. Our analyses and results using the concepts of process and state variables and their linear and nonlinear connectivity over TropSA shed new light on the fundamental role of soil moisture within the LAFs and the water and energy budgets over the region at interannual timescales.

Acknowledgments:
We would like to thank Kevin Trenberth (NCAR) for his valuable comments and suggestions of an earlier version of the manuscript.

Conflicts of Interest:
Authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.