Functional convergence of biosphere–atmosphere interactions in response to meteorological conditions

Understanding the dependencies of the terrestrial carbon and water cycle with meteorological conditions is a prerequisite to anticipate their behaviour under climate change conditions. However, terrestrial ecosystems and the atmosphere interact via a multitude of variables across temporal and spatial scales. Additionally these interactions might differ among vegetation types or climatic regions. Today, novel algorithms aim to disentangle the causal structure behind such interactions from empirical data. The estimated causal structures can be interpreted as networks, where nodes represent relevant meteorological 5 variables or land-surface fluxes, and the links the dependencies among them (possibly including time-lags and link strength). Here we derived causal networks for different seasons at 119 eddy-covariance flux tower observations in the FLUXNET network. We show that the networks of biosphere–atmosphere interactions are strongly shaped by meteorological conditions. For example, we find that temperate and high latitude ecosystems during peak productivity exhibit very similar biosphere– atmosphere interaction networks as tropical forests. In times of anomalous conditions like droughts though, both ecosystems 10 behave more like typical Mediterranean ecosystems during their dry season. Our results demonstrate that ecosystems from different climate zones or vegetation types have similar biosphere–atmosphere interactions if their meteorological conditions are similar. We anticipate our analysis to foster the use of network approaches as they allow a more comprehensive understanding of the state of ecosystem functioning. Long term or even irreversible changes in network structure are rare and thus can be indicators of fundamental functional ecosystem shifts. 15

correlation between two variables X and Y given a variable set Z is defined as the correlation between the residuals of X and Y after regressing out the (potentially multivariate) conditions Z. The conditions Z can consist of lagged third variables or time-lags of X and Y .

90
PCMCI has two phases. In the first phase, the 'condition selection', a superset of lagged parents (up to some maximum time lag τ max ) of each variable, X j t , is estimated based on a fast variant of the PC algorithm (Spirtes and Glymour, 1991). A parent of X j t is any lagged variable, X i t−τ , that is directly influencing X j t . This can be the own past, i = j, τ > 0 or other variables, i = j, τ > 0. A pseudo-code of this procedure is given in the supplementary materials of Runge et al. (2019b). In the second phase, 'momentary conditional independence' (MCI) is estimated among all pairs of contemporaneous and lagged variables 95 (X i t−τ , X j t ) for τ ≥ 0. The MCI test removes the influence of the lagged drivers (obtained in the first phase) using ParCorr and yields link strengths and p-values (based on a two-sided t-test). The link strength is here given by the MCI partial correlation .
In short, the MCI value gives an estimate of dependence between two time series, one potentially lagged, with the influence of other lagged drivers including autocorrelation removed, yielding a better interpretation of the strength of a causal mechanism than the common Pearson correlation. For a more detailed discussion of the interpretation, see Runge et al. (2019b). As a 100 particular partial correlation, the MCI value is independent of the variables' mean value and is normalised in [-1, 1] and can, hence, be compared between variable pairs with different units of measurement. Lagged links are directed forward in time.
Contemporaneous dependencies are left undirected as no time information reveals the direction of influence unless they are defined as unidirectional by the user (pcmci parameter selected_links, see table B1). A causal interpretation of links rests on the standard assumptions of causal discovery. Here we assume time order, the causal Markov condition, faithfulness, causal 105 sufficiency, causal stationarity, and no contemporaneous causal effects. The use of ParCorr additionally requires stationarity in the mean and variance and linear dependencies (Runge et al., 2019b). In particular, a statistical independence (here at a 0.1 two-sided significance level) between a pair of variables conditional on the other lagged variables is interpreted as the absence of a causal link (Faithfulness condition). On the other hand, a causal interpretation of the estimated links is here to be understood only with respect to the variables included in the analysis. The dependence structure among variables can finally 110 be visualised by weighted networks with the nodes representing the variables and the links significant dependencies with its strengths given by the MCI partial correlation.

Network Estimation
Dependencies are estimated using PCMCI among the variables R g , T, NEE, VPD, H and LE using time lags ranging from zero to five. As was already discussed by Krich et al. (2020), eddy covariance data and the choice of our variable set do not links while leaving the detection of true links largely unaffected. We estimated networks in sliding windows of three months, taking the centre month as the time index of each network. The sliding windows help to capture the temporal evolution of biosphere-atmosphere interactions and provide enough data points for the network estimation via PCMCI. Additionally, we improve stationarity of the data further and address the requirement of causal stationarity, i.e. a causal link persists throughout 125 the time period of network estimation. Further we set R g as a potential driver of the system (by excluding its parents from the PCMCI parameter 'selected_links', see table B1). We acknowledge the possibility of R g being influenced by other variables, e.g. via transpiration and subsequent cloud formation. Yet, on the ecosystem scale we work with, we presume this effect to be rather small and likely dominated by lateral transport. Besides these possibilities, setting R g as driver can account for remaining non stationarities (Runge, 2018). Missing data was flagged as such and is ignored by PCMCI. To avoid effects on the network 130 structure from gap-filling we used the following quality flag thresholds. A daily datapoint is not used if its quality flag is below 0.6 (i.e. more than 60% of measured and good quality gap filled data). In case more than 25% of datapoints of the three month window are flagged as bad quality, the time window is removed from the analysis. To analyse the factors influencing network structure, we consider the mean values over the respective time period of the variables included in the network calculation, and additionally those of GPP, P and SWC. GPP, P and SWC were not included in the network calculation because certain 135 characteristics can impinge on network estimation. GPP is derived using NEE and T. Any of the links GPP-T and GPP-NEE thus could be due to its processing rather than an actual dependence. P, on the other hand, typically yields non intuitive results due to its binary character (precipitation of a certain amount -zero precipitation), while its effects occur more smoothly (e.g. increase in transpiration or respiration) and its strong deviation from a normal distribution. Further, it can happen that over the time period of network estimation no precipitation occurs rendering such periods not analysable. The issue with SWC is 140 its lower availability and for those sites that have such measurements it might be applied at differing depth. The depth that is mostly present is at shallow depths of 5 or 10 cm. The upper soil layer, however, dries out quickly and can explain only little of the latent heat flux.

Dimensionality Reduction
For the dimensionality reduction, we tested principal component analysis (PCA; Pearson, 1901), t-distributed stochastic neigh-145 bour embedding (t-SNE; Maaten and Hinton, 2008), and uniform manifold approximation and projection (UMAP;McInnes et al., 2018). PCA is the standard method for dimensionality reduction, it is commonly used, linear, fast, and easily interpretable regarding the meaning of its axes (the principal components). A PCA embedding typically fails to reveal complex clusterings, because it maintains large scale gradients but often produces embeddings in which far away points appear very close in the embedding. In contrast t-SNE aims to preserve local neighbourhoods. Therefore it calculates first similarity scores 150 for each point pair using euclidean distances and Gaussian distributions. Subsequently it randomly projects the data onto the lower dimensional space and attempts to rearrange points in a way that the previously determined similarities are obtained.
To assess the similarities in the low dimensional space, however, it uses a Student-t distribution. This helps to separate points which are also originally separated. This procedure makes t-SNE very good at visualising clusters in the data and non-linear relationships. Drawbacks are the difficult interpretability of the embedding axes due to the non-linear nature and its fairly long 155 computation time for large datasets. Further, distances between far separated points and those belonging to different clusters in the embedding space are not (necessarily) comparable to the original distances. This is as t-SNE does not preserve both the global and local structure at the same time, which is attempted by UMAP. UMAP was developed as an improvement of t-SNE regarding structure preservation and results also in a shorter run time especially for higher dimensions. A comparison of t-SNE and UMAP is given in appendix C in McInnes et al. (2018). According to Kobak and Linderman (2019), the global structure 160 preservation of UMAP is not an inherent characteristic of the method itself but rather stems from the choosen initialization.
As we are dealing with an unsupervised method there is no obvious measure to assess the quality of an embedding, as each method optimises a different error function. A measure commonly used for the comparison and characterisation of dimensionality methods is the agreement between K-ary neighborhoods (the K nearest points to an observation) in the high dimensional and low dimensional space. The measure R NX (K) (Lee et al., 2015) gives a measure of the improvement of the embedding of 165 K-ary neighborhoods over random embeddings. For an embedding with random coordinates we obtain R NX (K) ≈ 0 and if the K-ary neighborhoods are perfectly preserved we obtain R NX (K) = 1. As this measure depends on the neighborhood size, K, we can draw a curve over K that characterizes if the method is better at maintaining global or local neighborhoods. The area under the R NX (K) curve gives an idea of the overall quality of the embedding. An intercomparison of the three dimensionality reduction methods using this measure shows t-SNE to perform best (see Fig. A1, B1, C1).

Distance Correlation
Distance correlation (Székely et al., 2007) is a non linear measure to quantify the dependence between two vectors. It has been used successfully to assess the influence of variables on the low dimensional embedding (Kraemer et al., 2020b). Székely et al. (2007) details its empirical definition for a sample (X, Y) = {(X k , Y k ) : k = 1, ..., n} with X ∈ R p and Y ∈ R q as follows: where V 2 n (X, Y) is the empirical distance covariance with V 2 n (X, Y) = 1 n 2 n k,l=1 A kl B kl . A kl and B kl are distance matrices defined by Distance correlation can be used to quantify the dependence between two sets of observations of differing dimensionality. In our case these two vectors are firstly a link strength or a underlying quantity of the networks (1d) and secondly the networks' 180 position in the low dimensional embedding (2d). The resulting dependence value is used to rank the quantities in their ability to describe the structure of the low dimensional embedding.

Clustering and median network trajectories
On the reduced space we applied a clustering method named Ordering Points To Identify the Clustering Structure (OPTICS; Ankerst et al., 1999 which consisted of at least three datapoints (networks within the t-SNE space belonging to the same ecosystem, representing the same month, in at least three years). The monthly median network is the average of the networks lying on (≥ 3 networks) or in the inner hull (< 3 networks).

Workflow
Our restrictions on the data length and quality resulted in a selection of 119 FLUXNET sites (Fig. 1a). Applying above 200 described procedure we obtained 10.038 networks for the different months and sites. An example network estimated by PCMCI is shown in Fig. 1c. The strongest and most consistent links are contemporaneous, indicating that interactions happen on time scales shorter than the time resolution. While lagged common drivers are excluded, contemporaneous links can still be spurious due to contemporaneous confounding (see Sect. 2.2). Nevertheless, we focus our analysis on these 15 links, as they contain most information. This is done by performing the dimensionality reduction on contemporaneous links and neglecting the lagged 205 ones. The rational of employing a dimensionality reduction is the following. Each of the estimated networks constitutes one observation in a high dimensional space with a network's links spanning its axes (Fig. 1d). Projecting this high dimensional space onto two dimensions ( Fig. 1e)   and most persistent links are contemporaneous (i.e. undirected). Thus we limit our analysis to those links. d) Each three-month network can be interpreted as an observation in a 15-dimensional space (each contemporaneous link is one dimension). e) Dimensionality reduction projects all interaction networks into a two dimensional space preserving its local neighbourhood structure. Here any subsequent analysis and interpretation will be realised.
as secondly static values like climate class, vegetation type or location. The secondary quantities are used to find covariates of the low dimensional embedding that can help to explain its structure. In a next step, we cluster the low dimensional embedding to further understand to which network structures the gradients of link strength lead (see Sect. 2.4) and calculate the cluster's average networks (a simple mean). Up to this point (Sect. 3.1 and 3.2), we analysed the manifold of biosphere-atmosphere interactions and can address the first part of our hypothesis. As each point of the low dimensional embedding represents the biosphere-atmosphere interactions of a specific ecosystem at a specific time we can investigate the behaviour of specific 220 ecosystems (see Sect. 2.7). Therefore we look at the monthly median and annual trajectories of certain ecosystems (Sect. 3.3 and 3.4). This leads to the answer of the second part of our hypothesis.

Results and Discussions
3.1 Two-dimensional embedding of biosphere-atmosphere networks To find the most suitable dimensionality reduction method, we evaluated three different methods (PCA, t-SNE and UMAP) with 225 respect to their ability to project the high dimensional network space onto two dimensions. To compare the low-dimensional embedding spaces, we used the R NX (K) measure (see Sect. 2.4) which quantifies how well neighbourhoods are preserved when projecting the high dimensional space onto fewer dimensions. We found that t-SNE achieved the best projection, by best preserving both local and distant neighbourhoods (cf. Sect. 2.4, Fig. A1, B1). This is unexpected as UMAP is said to intentionally preserve the global structure. Yet, as can be seen in Fig. 4a, the networks almost form a continuum. Thus, by 230 maintaining the local neighbourhood structure, also the global structure is preserved within t-SNE.
The two-dimensional embedding by t-SNE of biosphere-atmosphere interactions is ordered primarily by dependencies including carbon flux (NEE) and energy distributions (LE, H). This can be seen in Fig. 2 which shows the 2d embedding colourcoded by the strength of individual links, i.e. MCI partial correlation values. The colouring reveals that the link strengths are ordered along gradients, i.e. they exhibit some dependence with the t-SNE axes. Using distance correlation to rank those gra-   as they were used as features in the dimensionality reduction, gradients in mean climatic and biospheric conditions were not.
This information thus must be entailed in the networks' structure. To better grasp the distribution of network structures, we further analyse the emerging clusters.

Clusters of characteristic ecosystem-atmosphere networks
As we apply a significance threshold to each link of the estimated network structures (see Sect. 2.3), the networks typically lack weak links. This leads to a certain degree of clustering among the networks, which we identified using the OPTICS approach (see Sect. 2.6; Ankerst et al., 1999)) ( Fig. 4a). Cluster boundaries are shown by the convex hulls in Fig. 4b, where we also visualise the mean networks of each cluster. This visualisation reveals that the mean networks of the clusters situated 255 at the embedding's edges can be regarded as archetypes of network structures, i.e. extremal, characteristic states (similar to the concept of endmember states). The four states can be described as follows: Type 1 is a sparsely connected network. Links, if present, are very weak and predominantly exist among atmospheric variables.
Mean atmospheric conditions are characterised by low energy input (low R g and T). Carbon and water fluxes are con-sequently close to zero, and daily averages of sensible heat can even reach negative values. Such conditions reflect high latitude ecosystem winter states experienced by ecosystems like the evergreen needle leaf forests (ENF) of Finnland, i.e. Hyytiälä (FI-Hyy) and Sodankyla (FI-Sod) as well as Canada, i.e., the UCI-1850 burn site (CA-NS1) and Quebec -Eastern Boreal (CA-Qcu) during December and January.
Type 2 consists of strong links among atmospheric variables but LE and NEE are weakly, not, or even negatively connected to the atmosphere, i.e. the meteorological variables. This network structure coincides with high energy input (high R g and 265 T) but low water availability (low P and SWC, high VPD). A high Bowen ratio, i.e. the ratio between sensible heat and latent heat, representing aridity, and low absolute carbon fluxes (

280
The archetypes of networks are located at the edges of the two-dimensional space and thus could define two imaginary axes.
From a physical point of view, energy is required for each process and interaction to occur, e.g. photosynthesis or evaporation (Bonan, 2015). Therefore, transitions along the axis connecting the network types 1 and 4 might be interpreted as energy controlled as dependencies among all variables fade or increase consistently. Transitions along the axis connecting network types 2 and 3 are explainable by a combination of water availability and a temperature gradient. Low water availability but 285 high temperatures cause shut down of stomatal conductance or ecosystems to enter a dormant state which leads to low carbon and water fluxes and low connectivity. On the other hand, sufficient water and medium temperatures (around the optimum of photosynthesis) allow for carbon and water fluxes but reduce the influence of varying temperatures leading to connected NEE and LE but disconnected T. And indeed these patterns and gradients exist. Mean R g is lowest at network type 1 and almost linearly increases towards network type 4. P is lowest at network type 1 and 2. In combination with high energy input network  increase from network type 1 to 4 (as radiation) but also from network type 3 to 2 and are actually rather low (8 • C to 15 • C) at network type 3 (see Fig 3). As meteorological conditions affect biosphere productivity, network type 1 and 2 exhibit low, type 3 medium and type 4 high productivity i.e. estimated as GPP. In short, the clustering revealed that changes in energy and 295 water availability can explain major transitions between different states of biosphere-atmosphere interactions. This is in line with a recent study showing that a variety of land-surface processes can be largely summarised by on the one hand productivity measures and on the other hand water and energy availability. Both, water and energy availability, need to be high for high productive states, yet the lack of either of them leads to low productivity (Kraemer et al., 2020a). This biosphere state triangle is found in our analysis by the network type 1 (cold -low connectivity), 2 (dry -NEE/LE weakly connected) and 4 (high 300 productivity -fully connected). Yet, a fourth network type (type 3) is naturally occurring in the t-SNE space as we here include interactions with the atmosphere.
Up to this point we have found strong evidence supporting our first hypothesis. The manifold of biosphere-atmosphere interactions can be represented rather well by two dimensions which we identified to be most consistent with energy and water availabilities. It is confined by four characteristic states and populated homogeneously by the observed network states. Having 305 an understanding of the low dimensional embedding's structure now allows us to analyse specific ecosystem behaviour.

Ecosystems' median trajectories
Each point in the reduced t-SNE space represents a biosphere-atmosphere interaction network for a given month and ecosystem. Hence, we can trace an ecosystem's trajectory through time. We are first focusing on an ecosystem's median monthly trajectory (see Sect. 2.7) within the low dimensional space. We can see that the median trajectories reflect seasonal patterns 310 of meteorological conditions (Fig. 5). For example, mid-latitude sites like FR-Pue (EBF), DE-Hai (DBF) and FI-Hyy (ENF) exhibit a strong seasonal variation of R g and span a long distance in the t-SNE space. In contrast, tropical ecosystems like GF-Guy (EBF) constantly have high R g and exhibit predominantly network type 4 indicative of high productive conditions -while DE-Hai or FI-Hyy reach this connectivity pattern only during peak growing season. US-SRM (WSA), however, has similar or even higher R g values throughout the year but barely manages to deviate from type 2 which is in accordance with 315 its low water availability. The amount of precipitation further aligns with differences and characteristics of the trajectories of FR-Pue, DE-Hai and FI-Hyy. For example, FI-Hyy shows some deviation towards edge 2 in February and March, FR-Pue in June, July and August. For both, mean precipitation is lowest during these months. These behaviours demonstrate what the previous figures ( Fig. 3 and 4) have already suggested: Ecosystem's populate the low dimensional space and migrate within as allowed by their climatic conditions. Thereby they can exhibit a wide range of interaction structures as can be seen from the 320 mid-latitude sites. As these behaviours are multi year averages they could resemble more ecosystem adaptation to median climatic conditions than flexible adjustment of biosphere-atmosphere interactions to quickly changing meteorological conditions.
If biosphere-atmosphere interactions are confined by adaptation shall be investigated in the final analysis section.

Deviations from ecosystem median trajectories
The remaining open question is, how flexible do the networks adjust to deviations from mean climatic conditions. Therefore,

325
we look at climatic anomalies. Figure 6 shows the trajectories of ecosystems during anomalous dry or wet conditions. During the European heatwave of 2003, in July and August the trajectories of two temperate central European forests, DE-Hai and DE-Tha, no longer manage to establish a network structure resembling network type 4, typical for these ecosystems during their high productive phase. Instead they are shifted towards network type 2, associated with drier conditions (Fig. 6a, b).
Similarly, the ecosystem BR-Sa3 (EBF) in the Brazilian tropical rainforest shows substantial deviations towards network type  (Fig. 6c). In contrast, US-Wkg is a grassland accustomed to dry conditions and thus predominantly exhibits low water and carbon fluxes resulting in network structures as of network type 2, i.e. water and carbon fluxes are barely or even disconnected. Carbon and water fluxes of semi-arid ecosystems, however, are known to respond quickly and strongly to sufficient precipitation (Potts et al., 2019;Leon et al., 2014;Reynolds et al., 2004). This sensitivity is found to carry over to the network structure as well. The network structure of US-Wkg becomes 335 fully connected (network type 4) in September 2014 with above average precipitation (NOAA) (Fig. 6d) (right). In winter month the Bowen ratio can turn negative. Nevertheless we set the lower limit of the y-axis to 0. As networks are calculated using a centred three month moving window, each month is ascribed a network. Thus, the behaviour of an ecosystem can be tracked by its monthly networks, which form trajectories for each year. An ecosystem's monthly median trajectory is composed of the two dimensional monthly median networks (see Sect. 2.7 for details).
range. The opposite seems to be valid. Biosphere-atmosphere interactions can follow flexibly atmospheric conditions and are not confined to certain states.

Functional convergence of biosphere-atmosphere interactions
We have seen that networks representing biosphere-atmosphere interactions strongly align with prevailing mean meteorological conditions. Moreover, the visualisation of ecosystem trajectories within the t-SNE space (Fig. 5, 6)  larger part of the biosphere-atmosphere interaction network indeed is a pure atmospheric network, i.e. Rg, T, VPD and H.
Thus strong associations of networks and their trajectories with atmospheric conditions could be dominated by changes in this atmospheric network. Fig. 2, however, suggests the opposite. The strongest gradients are given by the links NEE-LE and Rg-LE and transitions along the axis connecting type 2 and 3 (cf. Fig. 4) are dominated by changes in biosphere connectivity, i.e. LE and NEE.

355
In fact, the dominance of climatic drivers in controlling the temporal evolution of ecosystem functioning emerges also in other studies (Musavi et al., 2017;Schwalm et al., 2017;Kraemer et al., 2020a) as they showed that carbon fluxes are primarily controlled by climatic factors. Yet, these and others also show the role of biotic factors in shaping the responses of ecosystem processes to climatic variability. For example, Musavi et al. (2017) revealed in a global ecosystem study that species diversity and ecosystem age decrease inter annual variability of GPP. Similarly, Wagg et al. (2017) showed biodiversity to increase long-360 term stability of ecosystem productivity. In regional studies Wales et al. (2020) found the stability of net primary production to be affected by the kind and severity of disturbances. Tamrakar et al. (2018) showed that seasonal carbon fluxes were more sensitive to environmental conditions in a homogeneous forest compared to a heterogeneous one. It would be of interest to investigate, to which degree the effects of biotic factors also translates to the sensitivity of the network structure.
Furthermore, extreme heat and drought events (Sippel et al., 2018) or compound events in general (Zscheischler et al., 2020) 365 can severely disrupt ecosystem functions. The time of recovery from such disturbances is a crucial parameter in assessing ecosystem resilience. Schwalm et al. (2017) showed that the recovery time measured as the recovery in GPP is primarily influenced by climate but secondarily by biodiversity and CO 2 fertilisation. Assessing the recovery time via GPP already puts the ecosystem functioning into focus. The here presented framework, i.e. the sensitivity of an ecosystem's network structure to meteorological conditions, might be a valuable asset to study reaction time and strength to and recovery from extreme events as 370 it not only utilises one variable but the interactions of a set of variables, thereby capturing more comprehensively an ecosystem state. A drawback is the reduced temporal resolution (a certain time period of daily or even half hourly measurements is aggregated to one network) which can be offset by the here used moving window approach to a certain degree. Especially with regard to climatic extreme conditions in recent years with observed vegetation dieback in, for example, DE-Hai (Schuldt et al., 2020), further studies could also shed light on the role of adaptation in shaping biosphere-atmosphere interactions. Our study 375 suggest that adaptation to a lesser degree limits the range of possible interactions but enables to sustain and persist certain conditions for longer periods. The focus of further studies thus could be to elucidate the role of biotic factors in influencing ecosystem trajectories as well as the role of adaptation and the response to extreme events.

Limitations of the study
Finally, we would like to take a critical view on our analysis approach. As stated in Sect. 2.2, PCMCI might fail to identify 380 some spurious links due to the occurrence of contemporaneous confounders. Thus networks can not be interpreted causally but this does not severely hinder their value for the current analysis. In addition we include a rather limited set of variables into the network estimation. Thus we cannot and do not claim that ecosystems become fully alike under similar meteorological conditions. Yet, on the timescale investigated the data shows, that the interactions among the chosen set of variables can be described by very similar structures. Follow up studies might search for and include further biosphere variables. Currently, an 385 analysis of the biotic effects on the network structure is hampered because the t-SNE space is not metric. Thus, for instance, the effect of a drought with similar magnitude in a boreal and temperate forest cannot simply be compared by the deviation from their median trajectory.

Conclusions
We analysed the functional behaviour of a variety of ecosystems using the FLUXNET database of carbon, water, and energy 390 flux measurements. In particular, we examined the interaction structure between biosphere-atmosphere fluxes as well as atmospheric state variables using PCMCI, a method to estimate causal relationships from empirical time series under certain assumptions. Using non-linear dimensionality reduction, we find evidence supporting our hypothesis that the manifold of existing states is bound by few, i.e. four, archetypes of network states. They are characterised on the one hand by a fully connected and almost unconnected network structure and on the other hand by an antagonistic coupling of carbon and water flux with 395 temperature -when one is strongly coupled, the other is decoupled. The transitions between these states correlate well with gradients of meteorological drivers, i.e. radiation and water availability. The movement of an ecosystem within that space therefore strongly aligns with changes in meteorological conditions. This, however, also leads to similar behaviour under similar conditions for strongly contrasting ecosystems. For example, forests of mid or even high latitudes exhibit similar interaction structure as tropical forests given high radiation and water availability during summer. Yet, this state can also be reached by 400 predominantly dry ecosystems like steppe grasslands given sufficient precipitation. In contrast if productive ecosystems are struck by a severe drought, like central European ecosystems in 2003, the behaviour can adapt more to that of a Mediterranean ecosystem. Thus the second part of our hypothesis must be rejected. The analysis shows that the biosphere-atmosphere interaction structure can adapt flexibly to prevailing conditions and is widely independent of vegetation type and climatic region.
Such behaviour is strong evidence for functional convergence of ecosystems, i.e. their behaviour is determined by a low num-405 ber of key processes. For further studies, we suggest, to focus on the role of biotic factors as, for example, plant functional types, ecosystem age and adaptation. These factors could play crucial roles in understanding the ecosystem copying strategies to climatic extremes.