Journal cover Journal topic
Biogeosciences An interactive open-access journal of the European Geosciences Union
Journal topic
Short summary
To closely monitor the state of our planet, we require systems that can monitor the observation of many different properties at the same time. We create indicators that resemble the behavior of many different simultaneous observations. We apply the method to create indicators representing the Earth's biosphere. The indicators show a productivity gradient and a water gradient. The resulting indicators can detect a large number of changes and extremes in the Earth system.
Final-revised paper
BG | Articles | Volume 17, issue 9
Biogeosciences, 17, 2397–2424, 2020
Biogeosciences, 17, 2397–2424, 2020

Research article 05 May 2020

Research article | 05 May 2020

Summarizing the state of the terrestrial biosphere in few dimensions

Summarizing the state of the terrestrial biosphere in few dimensions
Guido Kraemer1,2,3,4, Gustau Camps-Valls2, Markus Reichstein1,3, and Miguel D. Mahecha1,3,4 Guido Kraemer et al.
  • 1Max Planck Institute for Biogeochemistry, Department for Biogeochemical Integration, 07745 Jena, Germany
  • 2Image Processing Laboratory, Universitat de València, 46980 Paterna (València), Spain
  • 3German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, 04103 Leipzig, Germany
  • 4Remote Sensing Centre for Earth System Research, Leipzig University, 04103 Leipzig, Germany

Correspondence: Guido Kraemer (


In times of global change, we must closely monitor the state of the planet in order to understand the full complexity of these changes. In fact, each of the Earth's subsystems – i.e., the biosphere, atmosphere, hydrosphere, and cryosphere – can be analyzed from a multitude of data streams. However, since it is very hard to jointly interpret multiple monitoring data streams in parallel, one often aims for some summarizing indicator. Climate indices, for example, summarize the state of atmospheric circulation in a region. Although such approaches are also used in other fields of science, they are rarely used to describe land surface dynamics. Here, we propose a robust method to create global indicators for the terrestrial biosphere using principal component analysis based on a high-dimensional set of relevant global data streams. The concept was tested using 12 explanatory variables representing the biophysical state of ecosystems and land–atmosphere fluxes of water, energy, and carbon fluxes. We find that three indicators account for 82 % of the variance of the selected biosphere variables in space and time across the globe. While the first indicator summarizes productivity patterns, the second indicator summarizes variables representing water and energy availability. The third indicator represents mostly changes in surface albedo. Anomalies in the indicators clearly identify extreme events, such as the Amazon droughts (2005 and 2010) and the Russian heat wave (2010). The anomalies also allow us to interpret the impacts of these events. The indicators can also be used to detect and quantify changes in seasonal dynamics. Here we report, for instance, increasing seasonal amplitudes of productivity in agricultural areas and arctic regions. We assume that this generic approach has great potential for the analysis of land surface dynamics from observational or model data.

1 Introduction

Today, humanity faces negative global impacts of land use and land cover change (Song et al.2018), global warming (IPCC2014), and associated losses of biodiversity (IPBES2019; Díaz et al.2019), to only mention the most prominent transformations. Over the past decades, new satellite missions (e.g., Berger et al.2012; Schimel and Schneider2019) along with the continuous collection of ground-based measurements (e.g., Wingate et al.2015; Nasahara and Nagai2015; Baldocchi2020) and the integration of both (Papale et al.2015; Babst et al.2017; Jung et al.2019) have increased our capacity to monitor the Earth's surface enormously. However, there are still large knowledge gaps limiting our capacity to monitor and understand the current transformations of the Earth system (Steffen et al.2015; Rosenfeld et al.2019; Yan et al.2019; Piao et al.2020).

Many recent changes due to increasing anthropogenic activity are manifested in long-term transformations. One prominent example is “global greening” that has been attributed to fertilization effects, temperature increases, and land use intensification (de Jong et al.2011; Zhu et al.2016; Piao et al.2019). It is also known that phenological patterns change in the wake of climate change (Schwartz1998; Parmesan2006). However, these phenological patterns vary regionally. In “cold” ecosystems one may find decreased seasonal amplitudes on primary production due to warmer winters (Stine et al.2009). Elsewhere, seasonal amplitude may increase in agricultural areas, for example, due to the so-called “green revolution” (Zeng et al.2014; Chen et al.2019). Another change in terrestrial land surface dynamics is induced by increasing frequencies and magnitudes of extreme events (Barriopedro et al.2011; Reichstein et al.2013). The consequences for land ecosystems have yet to be fully understood (Flach et al.2018; Sippel et al.2018) and require novel detection and attribution methods tailored to the problem (Flach et al.2017; Mahecha et al.2007a). While extreme events are typically only temporary deviations from a normal trajectory, ecosystems may change their qualitative state permanently, for example shift from grassland to shrubland. Such shifts or tipping points can be induced by changing environmental conditions or direct human influence, and they pose yet another problem that needs to be considered (Lenton et al.2008). The question we address here is how to uncover and summarize changes in land surface dynamics in a consistent framework. The idea is to simultaneously take advantage of a large array of global data streams, without addressing each observed phenomenon in a specific domain only. We seek to develop an integrated approach to uncover changes in the land surface dynamics based on a very generic approach.

The problem of identifying patterns of change in high-dimensional data streams is not new. Extracting the dominant features from high-dimensional observations is a well-known problem in many disciplines. One approach is to manually define indicators that are known to represent important properties such as the “Bowen ratio” (Bowen1926, find a more complete description of the concept in Sect. 3.3). Another one consists in using machine learning to extract unique, and ideally independent features from the data. In the climate sciences, for instance, it is common to summarize atmospheric states using empirical orthogonal functions (EOFs), also known as principal component analysis (PCA; Pearson1901). The rationale is that dimensionality reduction only retains the main data features, which makes them more easily accessible for analysis. One of the most prominent examples is the description of the El Niño–Southern Oscillation (ENSO) dynamics in the multivariate ENSO index (MEI; Wolter and Timlin2011), an indicator describing the state of the regional circulation patterns at a certain point in time. The MEI is a very successful index that can be easily interpreted and used in a variety of ways; most basically it provides a measure for the intensity and duration of the different quasi-cyclic ENSO events, but it can also be associated with its characteristic impacts, e.g., seasonal warming, changes in seasonal temperatures, and overall dryness in the Pacific Northwest of the United States (Abatzoglou et al.2014); drought-related fires in the Brazilian Amazon (Aragão et al.2018); and crop yield anomalies (Najafi et al.2019).

In plant ecology, indicators based on dimensionality reduction methods are used to describe changes to species assemblages along unknown gradients (Legendre and Legendre1998; Mahecha et al.2007a). The emerging gradients can be interpreted using additional environmental constraints, or based on internal plant community dynamics (van der Maaten et al.2012). It is also common to compress satellite-based Earth observations via dimensionality reduction to get a notion of the underlying dynamics of terrestrial ecosystems. For instance, Ivits et al. (2014) showed that one can understand the impacts of droughts and heat waves based on a compressed view of the relevant vegetation indices. In general, dimensionality reduction is the method of choice to compress high-dimensional observations in a few (ideally) independent components with little loss of information (Van Der Maaten et al.2009; Kraemer et al.2018).

Understanding changes in land–atmosphere interactions is a complex problem, as all aforementioned patterns of change may occur and interact: land cover change may alter biophysical properties of the land surface such as (surface) albedo with consequences for the energy balance (Song et al.2018). Long-term trends in temperature, water availability, or fertilization may impact productivity patterns and biogeochemical processes (Zhu et al.2016; Sitch et al.2015). In fact, these land surface dynamics have implications for multiple dimensions and require monitoring of biophysical state variables such as leaf area index, albedo, etc., as well as associated land–atmosphere fluxes of carbon, water, and energy.

Here, we aim to summarize these high-dimensional surface dynamics and make them accessible for subsequent interpretations and analyses such as mean seasonal cycles (MSCs), anomalies, trend analyses, breakpoint analyses, and the characterization of ecosystems. Specifically, we seek a set of uncorrelated, yet comprehensive, state indicators. We want to have a set of very few indicators that represent the most dominant features of the above-described temporal ecosystem dynamics. These indicators should also be uncorrelated, so that one can study the system state by looking and interpreting each indicator independently. The approach should also give an idea of the general complexity contained in the available data streams. If more than a single indicator is required to describe land surface dynamics accurately, then these indicators shall describe very different aspects. While one indicator may describe global patterns of change, others could be only relevant in certain regions, for certain types of ecosystems, or for specific types of impacts. The indicators shall have a number of desirable properties: (1) represent the overall state of observations comprising the system in space and time, (2) carry sufficient information to allow for reconstructing the original observations faithfully from these indicators, (3) be of much lower dimensionality than the number of observed variables, and (4) allow intuitive interpretations.

In this work, we first introduce a method to create such indicators, and then we apply the method to a global set of variables describing the biosphere. Finally, to prove the effectiveness of the method, we interpret the resulting set of indicators and explore the information contained in the indicators by analyzing them in different ways and relating them to well-known phenomena.

Muller et al. (2011)Martens et al. (2017)Martens et al. (2017)Disney et al. (2016)Tramontana et al. (2016)Jung et al. (2019)Tramontana et al. (2016)Jung et al. (2019)Tramontana et al. (2016)Jung et al. (2019)Martens et al. (2017)Tramontana et al. (2016)Jung et al. (2019)Martens et al. (2017)Tramontana et al. (2016)Jung et al. (2019)Muller et al. (2011)

Table 1Variables used describing the biosphere. For a description of the variables, see Appendix A.

Download Print Version | Download XLSX

2 Methods

2.1 Data

Table 1 gives an overview of the data streams used in this analysis (for a more detailed description see Appendix A). For an effective joint analysis of more than a single variable, the variables have to be harmonized and brought to a single grid in space and time. The Earth System Data Lab (ESDL;, last access: 23 April 2020; Mahecha et al.2020) curates a comprehensive set of data streams to describe multiple facets of the terrestrial biosphere and associated climate system. The data streams are harmonized as analysis-ready data on a common spatiotemporal grid (equirectangular grid 0.25 in space and 8 d in time, 2001–2011), forming a 4D hypercube, which we call a “data cube”. The ESDL not only curates Earth system data, but also comes with a toolbox to analyze these data efficiently. For this study, we chose all available variables in the ESDL v1.0 (the most recent version available at the time of analysis), divided the available variables into meteorological and biospheric variables and discarded the atmospheric variables. We also discarded variables with distributions that are badly suited for a linear PCA (e.g., burned area contains mostly zeros) and variables with too many missing values. The only dataset that was added post hoc was fAPAR, which represents an important aspect of vegetation which was not available in the data cube at the time of analysis (it is part of the most recent version of the data cube).

The datasets taken from Tramontana et al. (2016) and Jung et al. (2019) are derived from flux tower measurements (Baldocchi2020). The flux towers are not equally distributed in climate space; i.e., there are many flux towers in temperate areas but much fewer in tropic and arctic regions, which may lead to less accurate data in these regions. These datasets also exclude large arid areas such as the Sahara and Gobi deserts and parts of the Arabian Peninsula which may affect the resulting loadings of the PCA slightly.

In this study, each variable was normalized globally to zero mean and unit variance to account for the different units of the variables, i.e., transform the variables to have standard deviations from the mean as the common unit. Because the area of the pixel changes with latitude in the equirectangular coordinate system used by the ESDL, the pixels were weighted according to the represented surface area. Only spatiotemporal pixels without any missing values were considered in the calculation of the covariance matrix.

2.2 Dimensionality reduction with PCA

As a method for dimensionality reduction, we used a modified principal component analysis to summarize the information contained in the observed variables. PCA transforms the set of d centered and, in this case, standardized variables into a subset of p, 1pd, principal components (PCs). Each component is uncorrelated with the other components, while the first PCs explain the largest fraction of variance in the data.

The data streams consist of d=12 observed variables at the same time and location. Each observation is defined in a d-dimensional space, xi∈ℝd, and we define the dataset by collecting all samples in the matrix X=[x1||xn]Rd×n. The observations are repeated in space and time and lie on a grid of lat×long×time. In our case, we have n=|lat|×|long|×|time|=720×1440×506=524,620,800 observations, where || denotes the cardinality of the dimension. Note that the actual number of observations was lower, n=106,360,156, because we considered land points only and removed missing values.

The fundamental idea of PCA is to project the data to a space of lower dimensionality that preserves the covariance structure of the data. Hence, the fundament of a PCA is the computation of a covariance matrix, Q. When all variables are centered to global zero mean and normalized to unit variance, the covariance matrix can in principle be estimated as

(1) Q = 1 n - 1 XX T = 1 n - 1 i = 1 n x i x i T .

However, in our case the data cube lies on a regular 0.25 grid and estimating Q as above would lead to overestimating the influence of dynamics in relatively small pixels of high latitudes compared to lower latitudes where each data point represents a larger area. Hence, one needs a weighted approach to calculate the covariance matrix,

(2) Q = 1 w i = 1 n w i x i x i T ,

where wi=cos (lati) and lati is the latitude of observation i, w=i=1nwi is the total weight, and n is the total number of observations. Equation (2) has the additional property that it can be computed sequentially on very big datasets, such as our Earth System Data Cube, by a consecutively adding observations to an initial estimate.

Note that the actual calculation of the covariance matrix is even more complicated, because summing up many floating-point numbers one by one can lead to large inaccuracies due to precision issues of floating-point numbers and instabilities of the naive algorithm (Higham1993; the same holds for the implementations of the sum function in most software used for numerical computing). Here, we used the Julia package WeightedOnlineStats.jl (, repository:, last access: 23 April 2020) (implemented by the first author of this paper), which uses numerically stable algorithms for summation, higher-precision numbers, and a map-reduce scheme that further minimizes floating-point errors.

Based on this weighted and numerically stable covariance matrix, the PCA can be computed using an eigendecomposition of the covariance matrix,

(3) Q = V Λ V T R d × d .

In this case, the covariance matrix Q is equal to the correlation matrix because we standardized the variables to unit variance. Λ is a diagonal matrix with the eigenvalues, λ1,,λd, in the diagonal in decreasing order and VRd×d, the matrix with the corresponding eigenvectors in columns. V can project the new incoming input data xi (centered and standardized) onto the retained PCs,

(4) y i = V T x i R d ,

where yi is the projection of the observation xi onto the d PCs.

The canonical measure of the quality of a PCA is the fraction of explained variance by each component, σi2, calculated as

(5) σ i 2 = λ i i = 1 d λ i .

To get a more complete measure of the accuracy of the PCA, we used the “reconstruction error” in addition to the fraction of explained variance. PCA allows a simple projection of an observation onto the first p PCs and a consecutive reconstruction of the observations from this p-dimensional projection. This is achieved by

(6) Y p = V p T X R p × n and X p = V p Y p R d × n ,

where Yp is the projection onto the first p PCs, Vp the matrix with columns consisting of the eigenvectors belonging to the p largest eigenvalues, and Xp the observations reconstructed from the first p PCs.

The reconstruction error, ei, was calculated for every point, xi, in the space–time domain based on the reconstructions from the first p principal components:

(7) e i = V p V p T x i - x i R d .

As this error is explicit in space, time, and variable, it allows for disentangling the contribution of each of these domains to the total error. This can be achieved by estimating the (weighed) mean square error,

(8) MSE = 1 w i w i e i 2 .

This approach can give a better insight into the compositions of the error than a single global error estimate based on the eigenvalues.

2.3 Pixel-wise analyses of time series

The principal components estimated as described above are ideally low-dimensional representations of the land surface dynamics that require further interpretation. These components have temporal dynamics that need to be understood in detail. One crucial question is how the dynamics of a system of interest deviate from its expected behavior at some point in time. A classical approach is inspecting the “anomalies” of a time series, i.e., the deviation from the mean seasonal cycle at a certain day of year.

Another key description of such system dynamics are trends. We estimated trends of the indicators as well as of their seasonal amplitude using the Theil–Sen estimator. The advantage of the Theil–Sen estimator is its robustness to up to 29.3 % of outliers (Theil1950a, b, c; Sen1968), while ordinary least-squares regression is highly sensitive to such values. The calculation of the estimator consists simply in computing the median of the slopes spanned by all possible pairs of points,

(9) slope i j = z i - z j t i - t j ,

where zi is the value of the response variable at time step i and ti the time at time step i. In our experiments, we computed the slopes separately per pixel and principal component with time as the predictor and the value of the principal component as the response variable.

To test the slopes for significance, we used the Mann–Kendall statistics (Mann1945; Kendall1970) and adjusted the resulting p values with the Benjamini–Hochberg method to control for the false discovery rate (Benjamini and Hochberg1995). Slopes with an adjusted p<0.05 were deemed significant.

To identify disruptions in trajectories, breakpoint detection provides a good framework for analysis. For the estimation of breakpoints, the generalized fluctuation test framework (Kuan and Hornik1995) was used to test for the presence of breakpoints. The framework uses recursive residuals (Brown et al.1975) such that a breakpoint is identified when the mean of the recursive residuals deviates from zero. We used the implementation in Zeileis et al. (2002). For practical reasons, here we only focus on the largest breakpoint.

Figure 1Example polygons and their areas, Eq. (10); the arrows indicate the directionality. (a) Clockwise polygon with a negative area. (b) Counterclockwise polygon with a positive area. (c) Chaotic polygon with a very low area. (d) Polygon with a single intersection and both a clockwise and counterclockwise portion. The clockwise portion is slightly larger than the counterclockwise portion; therefore the area is slightly negative.


The analysis of a different type of dynamic considers bivariate relations. In the context of oscillating signals it is particularly instructive to quantify their degree of phase shift and direction – even if both signals are not linearily related. A “hysteresis” would be such a pattern describing how the pathways AB and BA between states A and B differ (Beisner et al.2003). We estimated hysteresis by calculating the area inside the polygon formed by the mean seasonal cycle of the combinations of two components.

(10) Area = 1 2 i = 1 n x i ( y i + 1 - y i - 1 ) ,

where n=46, the number of time steps in a year, and xi and yi are the mean seasonal cycle of two PCs at time step i. The polygon is circular; i.e., the indices wrap around the edges of the polygon so that x0=xn and xn+1=x1. This formula gives the actual area inside the polygon only if it is non-self-intersecting and the vertices run counterclockwise. If the vertices run clockwise, the area is negative. If the polygon is shaped like an 8, the clockwise and counterclockwise parts will cancel each other (partially) out. Trajectories that have larger amplitudes will also tend to have larger areas as illustrated in Fig. 1.

3 Results and discussion

In the following, we first briefly present and discuss the quality of the global dimensionality reduction (Sect. 3.1) and interpret the individual components from an ecological point of view (Sect. 3.2). We summarize the global dynamics that we uncovered in the low-dimensional space (Sect. 3.3). We characterize the contained seasonal dynamics (Sect. 3.4), including spatial patterns of hysteresis (Sect. 3.5). We then describe global anomalies of the identified trajectories (Sect. 3.6) and discuss the identified anomalies in depth based on local phenomena (Sect. 3.7). Finally, we present global trends and their breakpoints (Sect. 3.7).

Figure 2(a) Fraction of explained variance of the PCA by component. The knee at component three suggests that components four and higher do not contribute much to total variance. (b) Rotation matrix of the global PCA model (also called loadings, Eq. 4). The columns of the rotation matrix describe the linear combinations of the (centered and standardized) original variables that make up the principal components. PC1 is dominated by primary-productivity-related variables, PC2 by variables describing water availability, and PC3 by variables describing albedo. Values of the rotation matrix are clamped to the range [0.5, 0.5]; the actual range of the values is [0.73, 0.74] and [0.46, 0.54] for the first three components.


3.1 Quality of the PCA

Figure 2a shows the explained fraction of variance (Eq. 5) for the global PCA based on the entire data cube. The two leading components explain 73 % of the variance from the 12 variables; additional components contribute relatively little additional variance (PC3 contributes 9 % and all subsequent PCs less than 7 %) each. This results in a “knee” at component 3, which suggests that two indicators are sufficient to capture the major global dynamics of the terrestrial land surface, but we will also consider the third components in the following analyses (Cattell1966).

Figure 3Reconstruction error of the data cube using varying numbers of principal components aggregated by the mean squared error. Reconstruction errors aggregated over all time steps and variables are shown in the left column: (a) using only the first component, (c) using the first two, (e) and using the first three. Corresponding right plots (b, d, f) show the mean reconstruction error aggregated by latitude.

We estimated the reconstruction error sequentially up to the first three principal components (Fig. 3). Regions that do not fit the model well show a higher reconstruction error. Considering one component only, the highest reconstruction errors appear in high latitudes but decrease strongly with each additional component and nearly vanish if the third component is included.

3.2 Interpretation of the PCA

The first PC summarizes variables that are closely related to primary productivity (GPP, LE, NEE, fAPAR) and therefore are highly interrelated (see Fig. 2b). The energy for photosynthesis comes from solar radiation, and fAPAR is an indicator for the fraction of light used for photosynthesis. The available photosynthetic radiation is used by photosynthesis to fix CO2 and to produce sugars that maintain the metabolism of the plant. The total uptake of CO2 is reflected in GPP, which is also closely related to water consumption. The flow of water within the plant is not only essential to enable photosynthesis but also drives the transport of nutrients from the roots. The uplift of water in the plant is ultimately driven by transpiration – together with evaporation from soil surfaces one can observe the integrated latent energy needed for the phase transition (LE). However, ecosystems also respire; CO2 is produced by plants in energy-consuming processes as well as by the decomposition of dead organic materials via soil microbes and other heterotrophic organisms. This total respiration can be observed as terrestrial ecosystem respiration (TER). The difference between GPP and TER is the net ecosystem exchange (NEE) rate of CO2 between ecosystems and the atmosphere (Chapin et al.2006). GPP and TER are also well represented in the first dimension.

The second component represents variables related to the surface hydrology of ecosystems (see Fig. 2b). Surface moisture, evaporative stress, root-zone soil moisture, and sensible heat are all essential indicators for the state of plant-available water. While surface moisture is a rather direct measure, evaporative stress is a modeled quantity summarizing the level of plant stress: a value of zero means that there is no water available for transpiration, while a value of 1 means that transpiration equals the potential transpiration (Martens et al.2017). Root-zone soil moisture is the moisture content of the root zone in the soil, the moisture directly available for root uptake. If this quantity is below the wilting point, there is no water available for uptake by the plants. Sensible heat is the exchange of energy by a change in temperature; if there is enough water available, then most of the surface heat will be lost due to evaporation (latent heat), and with decreasing water availability more of the surface heat will be lost due to sensible heat, making this an indicator of dryness as well.

We observe that the third component is most strongly related to albedo (Fig. 2b). Albedo describes the overall reflectiveness of a surface. Here we refer to broadband (400–3000 nm) surface albedo; for an exact definition see Appendix A. Light surfaces, such as snow and sand, reflect most of the incoming radiation, while surfaces that have a high liquid water content or active vegetation absorb most of the incoming radiation. Local changes to albedo can be due to many causes, e.g., snowfall, vegetation greening and browning, or land use change.

The relation of PC3 to productivity and hydrology is opposite to what we would expect from an albedo axis. Because vegetation uses radiation as an energy source, albedo is negatively correlated with the productivity of vegetation, hence the negative correlation of albedo with PC1. Given that water also absorbs radiation, we can observe a negative correlation of albedo with PC2 (see Fig. 2b). We observe that PC1 and PC2 are positively correlated with PC3 on the positive portion of their axes (see Fig. 4d and f), which means counterintuitively that the index representing albedo is positively correlated with primary productivity and moisture content. Finally we can observe that PC1 and PC2 have a much higher reconstruction error in snow-covered regions, which is strongly improved by adding PC3 (see Fig. 3f). Therefore the third component should be regarded mostly as a binary variable that introduces snow cover, as the other information that is usually associated with albedo is already contained in the first two components.

Figure 4Trajectories of some points (colored lines) and the area-weighted density over principal components one and two (the gray background shading shows the density) for (a, c, e) the raw trajectories and (b, d, f) the mean seasonal cycle. The trajectories are shown in the space of PC1–PC2 (first row), PC1–PC3 (second row), and PC2–PC3 (third row). The trajectories were chosen to cover a large area in the space of the first two principal components. Some of the trajectories have an arrow indicating the direction. The numbers illustrate the value of some variables; for units see Table 1. Description of the points is as follows. Red: tropical rain forest, 2.625 S, 67.625 W; blue: maritime climate, 52.375 N, 7.375 E; green: monsoon climate, 22.375 N, 82.375 E; purple: subtropical, 34.875 N, 117.625 W; orange: continental climate, 52.375 N, 44.875 E; yellow: arctic climate, 72.375 N, 119.875 E.


3.3 Distribution of points in PCA space

The bivariate distribution of the first two principal components forms a “triangle” (gray background in Fig. 4a). At the high end of PC1 we find one point of the triangle in which ecosystems have a high primary productivity (high values of GPP, fAPAR, LE, TER, and evaporation), mostly limited by radiation. On the lower end of the first principal component we find the other two points of the triangle describing two alternative states of low productivity. These can happen either when the second principal component coincides with temperature limitation (the negative extreme of the second principal component) as seen in the lower left corner of the distribution in Fig. 4a and b or due to water limitation (positive extreme of the second principal component, the upper left corner in Fig. 4a). This pattern reflects the two essential global limitations of GPP in terrestrial ecosystems (Anav et al.2015).

Both components form a subspace in which most of the variability of ecosystems takes place. Component one describes productivity and component two the limiting factors to productivity. Therefore, we can see that most ecosystems with high values on component one (a high productivity) are at the approximate center of component two. When ecosystems are found outside the center of component two, they have lower values on component one (lower productivity) because they are limited by water or temperature (see Fig. 4b).

Figure 5The background shading shows the distribution of the mean seasonal cycle of the spatial points (see Fig. 4). The contour lines represent the reconstruction of the variables from the first two principal components. The reconstructed variables are (a) latent heat (LE), (b) sensible heat (H), and (c) log10SensibleHeatLatentHeat, the log 10 of the Bowen ratio. Note that the LE and H have been considered in the construction of the PCs and hence are a linear function of the PCs. The Bowen ratio, instead, was not considered here and clearly responds in a nonlinear form.


To further interpret the triangle we analyze how the Bowen ratio embeds in the space of the first two dimensions. Energy fluxes from the surface into the atmosphere can represent either a radiative transfer (sensible heat) or evaporation (latent heat). Their ratio is the “Bowen ratio”, B=HLE, (Bowen1926; see also Fig. 5). When water is available most of the available energy will be dissipated by evaporation, B<1, resulting in a high latent heat flux. Otherwise, the transfer by latent heat will be low and most of the incoming energy has to be dissipated via sensible heat, B>1. In higher latitudes, there is relatively limited incoming radiation and temperatures are low; therefore there is not much energy to be dissipated and both heat fluxes are low. A high sensible heat flux is an indicator of water limitation.

3.4 Seasonal dynamics

The leading principal components represent most of the variability of the space spanned by the observed variables, summarizing the state of a spatiotemporal pixel efficiently. This means that the PCs track the state of a local ecosystem over time (Fig. 4a) or, in the case of the mean seasonal cycle, time of the year (Fig. 4b). For a representation of the state of the first three components in time and space, see Appendix Fig. B1.

A first inspection reveals a substantial overlap of seasonal cycles of very different regions of the world. We also see that very different ecosystems may reach very similar states in the course of the season, even though their seasonal dynamics are very different. For instance, a midlatitude pixel (blue trajectory in Fig. 4) shows very similar characteristics to tropical forests during peak growing season. This indicates that an ecosystem of the midlatitudes can reach similar levels of productivity and water availability as a tropical rain forest (see also Appendix Fig. C1). Likewise, for the first two components, many high-latitude areas show similar characteristics to midlatitude areas during winter (low latent and sensible energy release as well as low GPP), and many dry areas such as deserts show similar characteristics to areas with a pronounced dry season, e.g. the Mediterranean.

Depending on their position on Earth, ecosystem states can shift from limitation to growth during the year (Fig. 4b, e.g. Forkel et al.2015). For example, the orange trajectory in Fig. 4, an area close to Moscow, shifts from a temperature-limited state in winter to a state of very high productivity during summer. Other ecosystems remain in a single limitation state with only slight shifts, such as the red trajectory in Fig. 4. In the corner of maximum productivity of the distribution, we find tropical forests characterized by a very low seasonality. We also observe that very different ecosystems can have very similar characteristics during their peak growing season; e.g. green (located in northeast India), blue (northwest Germany), and orange (located close to Moscow) trajectories have very similar characteristics during peak growing season compared to the red trajectory.

The third component shows a different picture. Due to a consistent winter snow cover in higher latitudes, the albedo is much higher and the amplitude of the mean seasonal cycle is much larger than in other ecosystems. Other areas show comparatively little variance on the third component and their relation to productivity and moisture content is even positively correlated to the third component, which is the opposite of what is expected from an albedo axis.

Figure 6Mean seasonal cycle of the first three principal components (in columns) during the seasons (in rows). Left column: first principal component. Middle column: second principal component. Right column: third principal component. Rows from top to bottom: equally spaced intervals during the year. Values have been clamped to 0.7 times their range to increase contrast.

The global pattern of the first principal component follows the productivity cycles during summer and winter (Fig. 6, left column) of the Northern Hemisphere, with positive values (high productivity, green) during summer and negative values (low productivity, brown) during winter. The tropics show high productivity all year. The global pattern shows the well-known green wave (Schwartz1994, 1998) because the first dimension integrates over all variables that correlate with plant productivity.

The second principal component (Fig. 6, middle column) tracks water deficiency: red and light red areas indicate water deficiency, light blue areas excess water, and dark blue areas water growth limitation due to cold. Areas which are temperature limited during winter but have a growing season during summer, such as boreal forests, change from dark blue in winter to light blue during the growing season. Areas which have low productivity during a dry season change their coloring from red to light red during the growing season, e.g the northwest of Mexico and southwest of the United States.

The third principal component (Fig. 6, right column) tracks surface reflectance. Therefore we can see the highest values in the arctic region during winter, and other areas vary much less in their reflectance throughout the year. Again, the third component shows a counterintuitive behavior in the midlatitudes, as it is positively correlated with productivity and therefore shows the opposite behavior of what would be expected from an indicator tracking albedo.

Although the principal components are globally uncorrelated, they covary locally (see Fig. D1). Ecosystems with a dry season have a negative covariance between PC1 and PC2, while ecosystems that cease productivity in winter have a positive covariance. Cold arid steppes and boreal climates show a negative covariance between PC1 and PC3. While other ecosystems that have a strong seasonal cycle show a positive correlation, many tropical ecosystems do not show a large covariance. A very similar picture is painted between the covariance of PC2 and PC3: boreal and steppe ecosystems show a negative covariance, while most other ecosystems show a more or less pronounced positive covariance, again depending on the strength of the seasonality.

Observing the mean seasonal cycle of the principal components gives us a tool to characterize ecosystems and may also serve as a basis for further analysis, such as a global comparison of ecosystems (Metzger et al.2013; Mahecha et al.2017).

Figure 7The area inside the mean seasonal cycles of (a) PC1–PC2, (b) PC1–PC3, and (c) PC2–PC3. The area is positive if the direction is counterclockwise and negative if the direction is clockwise. Most of the trajectories need a strong seasonal cycle to show a pronounced hysteresis effect. If the mean seasonal cycle intersects, the areas cancel each other out, e.g. the green trajectory of Fig. 4b.

3.5 Hysteresis

The alternative return path between ecosystem states forming the hysteresis loops arises from the ecosystem tracking seasonal changes in the environmental condition, e.g. summer–winter or dry–rainy seasons (Fig. 4b). Hysteresis is a common occurrence in ecological systems (Folke et al.2004; Blonder et al.2017; Renner et al.2019). For instance, a hysteresis loop can be found when plotting soil respiration against soil temperature (Tang et al.2005). The sensitivity of soil respiration to soil temperature changes seasonally due to changing soil moisture and photosynthesis (by supplying carbon to the rhizosphere), producing a seasonally changing hysteresis effect (Gaumont-Guay et al.2006; Richardson et al.2006; Zhang et al.2018). Biological variables also show a hysteresis effect in their relations with atmospheric variables; e.g. Mahecha et al. (2007b) found a hysteresis effect between seasonal NEE, temperature, and a number of other ecosystem and climate-related variables. Here we look at the mean seasonal cycles of pairs of indicators and the area they enclose.

The orange trajectory (area close to Moscow) in Fig. 4b shows that the paths between maximum and minimum productivity can be very different, in contrast to the blue trajectory located in the northwest of Germany which also has a very pronounced yearly cycle but shows no such effect. Figure 4 also indicates that the area inside the mean seasonal cycles of PC1–PC2 and PC1–PC3 shows important characteristics while hysteresis in PC2–PC3 is a much less pronounced feature; i.e., we can only see a pronounced area inside the yellow curve in Fig. 4f.

The trajectories that show a more pronounced counterclockwise hysteresis effect in PC1–PC2 (Fig. 7a) are areas with a warm and temperate climate and partially those that have a snow climate with warm summers, i.e., areas that have pronounced growing, dry, and wet seasons and therefore shift their limitations more strongly during the year. That means the moisture reserves are depleted during growing season, and therefore the return path has higher values on the second principal component (the climatic zones are taken from the Köppen–Geiger classification; Kottek et al.2006). We can also see that areas with dry winters tend to have a clockwise hysteresis effect, e.g. many areas in East Asia. Due to the humid summers there is no increasing water limitation during the summer months which causes a decrease for PC2 instead of an increase. Other areas with clockwise hysteresis can be found in winter dry areas in the Andes and the winter dry areas north and south of the African rain forests. Tropical rain forests do not show any hysteresis effect due to their low seasonality. In general we can say that the area inside the mean seasonal cycle trajectory of PC1–PC2 depends mostly on water availability in the growing and non-growing seasons, i.e., the contrast of wet summer and dry winter vs. dry summer and wet winter.

The hysteresis effect on PC1–PC3 (Fig. 7b) shows a pronounced counterclockwise MSC trajectory mostly in warm temperate climates with dry summers, while it shows a clockwise MSC trajectory in most other areas; again tropical rain forests are an exception due to their low seasonality. The most pronounced clockwise MSC trajectories can be found in tundra climates in arctic latitudes, where we have a consistent winter snow cover and a very short growing period. A counterclockwise rotation can be found in summer dry areas, such as the Mediterranean and California, but also some more humid areas, such as the southeast United States and the southeast coast of Australia. In these areas we can find a decrease for PC3 during the non-growing phase which probably corresponds to a drying out of the vegetation and soils.

The hysteresis effect on PC2–PC3 (Fig. 7c) mostly depends on latitude. There is a large counterclockwise effect in the very northern parts, due to the large amplitude of PC3. The amplitude gets smaller further south until the rotation reverses in winter dry areas at the northern and southern extremes of the tropics and disappears at the equatorial humid rain forests.

We can see that the hysteresis of pairs of indicators represents large-scale properties of climatic zones. The enclosed area and the direction of the rotation provide interesting information. Hysteresis can provide information on the seasonal availability of water, seasonal dry periods, or snowfall. With the method presented here, we can not observe intersecting trajectories, which would probably provide even more interesting insights (e.g. the green trajectory in Fig. 4b).

Figure 8Anomalies of the first three principal components. The brown–green contrast shows the anomalies on PC1, a relative low productivity or greening, respectively. The blue–red contrast shows the anomalies on PC2, a relative wetness or dryness, respectively. The brown–purple contrast shows the anomaly on PC3, a relative deviation in albedo. Panels (a), (e), and (i) are maps showing the anomalies of PC1–PC3, respectively, on 1 January 2001. Panels (b), (c), and (d) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (a). The effects of the floods on the Horn of Africa (2006) and the Russian heat wave (2010) are highlighted by circles. Panels (f), (g), and (h) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (e). Strong droughts in the Amazon during 2005 and 2010 can be observed as large red spots on the fringes of the Amazon basin (highlighted by circles). Panels (j), (k), and (l) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (i). A strong snowfall event affecting central and southern China is marked as circles.

3.6 Anomalies of the trajectories

The deviation of the trajectories from their mean seasonal cycle should reveal anomalies and extreme events. These anomalies have a directional component which makes them interpretable the same way the original PCs are. Therefore one can infer the state of the ecosystem during an anomaly. For instance the well-known Russian heat wave in summer 2010 (Flach et al.2018) appears in Fig. 8 as a dark brown spot in the southern part of the affected area, indicating lower productivity, and as a thin green line in the northern parts, indicating increased productivity. This confirms earlier reports in which only the southern agricultural ecosystems were negatively affected by the heat wave, while the northern predominantly forest ecosystems rather benefited from the heat wave in terms of primary productivity (Flach et al.2018).

Another example of an extreme event that we find in the PCs is the very wet November rainy season of 2006 in the Horn of Africa after a very dry rainy season in the previous year. This event was reported to bring heavy rainfall and flooding events which caused an emergency for the local population but also increased ecosystem productivity (Nicholson2014). The rainfall event appears as green and blue spots in Fig. 8b and c, preceded by the drought events which appear as red and brown spots.

Figure 8f and g also show the strong drought events in the Amazon, particularly the droughts of 2005 and 2010 (Doughty et al.2015; Feldpausch et al.2016) appear strongly north and south of the Amazon basin. The central Amazon basin does not show these strong events, because the observable response of the ecosystem was buffered due to the large water storage capacity in the central Amazon basin.

Another extreme event that can be seen is the extreme snow and cold event affecting central and south China in January 2008, causing the temporary displacement of 1.7 million people and economic losses of approximately USD 21 billion (Hao et al.2011). This event shows up clearly on PC2 and PC3 as cold and light anomalies, respectively (see Fig. 8k and f).

Figure 9Trajectories of the first two principal components for single pixels. (a) Deforestation increases the seasonal amplitude of the first two PCs (Brazilian rain forest, 9.5 S, 63.5 W). The red line shows the trajectory before 2003 and the blue line the trajectory 2003 and later. A strong increase in seasonal amplitude can be observed after 2003. (b) The heat wave is clearly visible in the trajectory (red, Russian heat wave, summer 2010, 56 N, 45.5 E). (c) Rainfall in the short rainy season (November–December) influences agricultural yield and can cause flooding (extreme flooding after drought, November 2006, 3 N, 45.5 E). (d) The European heat wave in summer 2003 was one of the strongest on record (France, 47.2 N, 3.8 E). The mean seasonal cycle of the trajectories is shown in purple.

3.7 Single trajectories

Observing single trajectories can give insight into past events that happened at a certain place, such as extreme events or permanent changes in ecosystems. The creation of trajectories is an old method used by ecologists, mostly on species assembly data of local communities, to observe how the composition changes over time (e.g. Legendre et al.1984; Ardisson et al.1990). In this context, we observe how the states of the ecosystems inside the grid cell shift over time, which comprises a much larger area than a local community but is probably also less sensitive to very localized impacts than a community-level analysis. One of the main differences of the method applied here from the classical ecological indicators is that the trajectories observed here are embedded into the space spanned by a single global PCA, and therefore we can compare a much broader range of ecosystems directly.

Figure 10(a, c, e) Trends in PC1–PC3, respectively (2001–2011). (b, d, f) Bivariate distribution of trends. Trends were calculated using the Theil–Sen estimator. Panels (a), (c), and (e) show significant trends only (p<0.05, Benjamini–Hochberg adjusted).

The seasonal amplitude of the trajectory in the Brazilian Amazon increases due to deforestation and crop growth cycles. Figure 9a shows an area in the Brazilian Amazon in Rondônia (9.5 S, 63.5 W) which was affected by large-scale land use change and deforestation. It can be seen that the seasonal amplitude increases strongly after the beginning of 2003. This increased amplitude could be due to any of the following reasons or a combination of them: deforestation decreases water storage capability and dries out soils, causing larger variability in ecosystem productivity. Therefore, during periods of no rain, large-scale deforestation can cause a shift in local-scale circulation patterns, causing lower local precipitation (Khanna et al.2017). Crop growth and harvest cause an increased amplitude in the cycle of productivity. An analysis of the trajectory can point to the nature of the change; however finding the exact causes for the change requires a deeper analysis.

The 2010 Russian heat wave has a very clear signal in the trajectories. Figure 9b shows the deviation of the trajectory during the Russian heat wave (red line) in an area east of Moscow (56 N, 45.5 E). In the southern grass- and croplands, the heat wave caused the productivity to drop significantly during summer due to a depletion of soil moisture. In the northern forested parts affected, the heat wave caused an increase in ecosystem productivity during spring due to higher temperatures combined with sufficient water availability. This shows the compound nature of this extreme event (see Fig. 8a and Flach et al.2018). The analysis of the trajectory points directly towards the different types of extremes and responses that happened in the biosphere during the heat wave.

Variability of rainfall during the November rainy season in the Horn of Africa (3 N, 45.5 E, Fig. 9c) shows the trajectory and points in November of the observed time. The November rain has implications for food security because the second crop season depends on it. In 2006, the rainfall events were unusually strong and caused widespread flooding and disaster but also higher ecosystem productivity (see also Fig. 8). This was especially devastating because it followed a long drought that caused crop failures. Note also the two rainy seasons in the mean seasonal cycle (purple line in Fig. 9c).

The 2003 European heat wave is reflected in the trajectories just like the 2010 Russian heat wave. Figure 9d shows the trajectory during the August 2003 heat wave in Europe (France, 47.2 N, 3.8 E). The heat wave was unprecedented and caused large-scale environmental, health, and economic losses (Ciais et al.2005; García-Herrera et al.2010; Miralles et al.2014). The 2010 heat wave was stronger than the 2003 heat wave but the strongest parts of the 2010 heat wave were in eastern Europe (see Fig. 8), while the center of the 2003 heat wave was located in France.

As we have seen here, observing single trajectories in reduced space can give us important insights into ecosystem states and changes that occur. While the trajectories can point us towards abnormal events, they can only be the starting points for deeper analysis to understand the details of such state changes.

3.8 Trends in trajectories

The accumulation of CO2 in the atmosphere should cause an increase in global productivity of plants due to CO2 fertilization, while larger and more frequent droughts and other extremes may counteract this trend. Satellite observations and models have shown that during the last decades the world's ecosystems have greened up during growing seasons. This is explained by CO2 fertilization, nitrogen deposition, climate change, and land cover change (Zhu et al.2016; Huang et al.2018; Anav et al.2015). Tropical forests especially showed strong greening trends during the growing season.

General patterns of trends that can be observed are a positive trend (higher productivity) on the first principal component in many arctic regions. Many of these regions also show a wetness trend, with the notable exception of the western parts of Alaska, which have become drier. This is important, because wildfires play a major role in these ecosystems (Jolly et al.2015; Foster et al.2019). These changes are also accompanied by a decrease for PC3 due to a loss in snow cover. A large-scale dryness trend can also be observed across large parts of western Russia. Increasing productivity can also be observed for large parts of the Indian subcontinent and eastern Australia. Negative trends in the first component can also be observed: they are generally smaller and appear in regions around the Amazon and the Congo Basin, but also in parts of western Australia. The main difference from previous analyses on the observations presented here is that Zhu et al. (2016), for example, looked only at trends during the growing season, while this analysis uses the entire time series to calculate the slope.

In the Amazon basin, we find a dryness trend accompanied by a decrease in productivity and a slight increase in PC3. In the Congo Basin, we find a wetness trend and an increasing productivity in the northern parts, while the southern part and woodland south of the Congo Basin show a strong dryness trend with decreased productivity. This is different to the findings of Zhou et al. (2014), who found a widespread browning of vegetation in the entire Congo Basin for the April–May–June seasons during the period 2000–2012. The findings of Zhou et al. (2014) are not reflected in our data, especially compared to the areas surrounding the Congo Basin. We can find only minor browning effects inside the basin, and our findings are more in line with the global greening (Zhu et al.2016), which shows a browning mostly outside the Congo Basin.

In eastern Australia we find a strong wetness and greenness trend which is due to Australia having a “millennium drought” since the mid-1990s with a peak in 2002 (Nicholls2004; Horridge et al.2005) and extreme floods in 2010–2011 (Hendon et al.2014).

Large parts of the Indian subcontinent show a trend towards higher productivity and an overall wetter climate. The greening trend in India happens mostly over irrigated cropland. However browning trends over natural vegetation have been observed but do not emerge in our analysis (Sarmah et al.2018). A very notable greening and wetness trend can be observed in Myanmar due to an increase in intense rainfall events and storms, although the central part experienced some strong droughts at the same time (Rao et al.2013). In Myanmar we also find one of the strongest trends in PC3 outside of the Arctic.

In large parts of the Arctic, a trend towards higher productivity can be observed. Vegetation models attribute this general increase in productivity to CO2 fertilization and climate change. The changes also cause changes to the characteristics of the seasonal cycles (Forkel et al.2016). Stine et al. (2009) found a decreased seasonal amplitude of surface temperature over northern latitudes due to winter warming.

The seasonal amplitude of atmospheric CO2 concentrations has been increasing due to climate change, causing longer growing seasons and changing vegetation cover in northern ecosystems (Forkel et al.2016; Graven et al.2013; Keeling et al.1996). Therefore we checked for trends in the seasonal amplitude, but because each time series only consists of 11 values (one amplitude per year), after adjusting the p values for false discovery rate, we could not find a significant slope. However, there were many significant slopes with the unadjusted p values; see the appendix, Fig. E1.

Another way to detect changes to the biosphere consists in the detection of breakpoints, which has been applied successfully to detect changes in global normalized difference vegetation index (NDVI) time series (de Jong et al.2011; Forkel et al.2013) or generally to detect changes in time series (Verbesselt et al.2010). A proof-of-concept analysis can be found in Fig. F1. We hope that applying this method to indicators instead of variables can detect a wider range of breakpoints analyzing a single time series.

3.9 Relations to other PCA-type analyses

One of the most popular applications of PCA in meteorology are EOFs, which typically apply PCA to a single variable, i.e., on a dataset with the dimensions lat×long×time, although EOFs can be calculated from multiple variables. EOFs can be calculated in S mode and R mode. If we matricize our data cube so that we have time in rows and lat×long×variables in columns, then S mode PCA works on the correlation matrix of the combined variable and space dimension. In T mode, the PCA works on the correlation matrix formed by the time dimension (Wilks2011). The PCA presented here works slightly differently. (1) We performed a different matricization (lat×long×time in rows and variables in columns) and then (2) the PCA works on the correlation matrix formed by the variables. Therefore in this framework we could call this a V mode PCA.

Ecological analyses usually use PCA with matrices of the shape object×descriptors. When calculating the PCA on the correlation matrix formed by the objects, then it is called a Q mode analysis. When the PCA is applied to the correlation matrix formed by the variables, then it is called an R mode analysis (Legendre and Legendre1998). The PCA carried out in this study is closest to an R mode analysis. In the present case the descriptors are the various data streams and the objects are the spatiotemporal pixels.

Using PCA as a method for dimensionality reduction means that we are assuming linear relations among features. A nonlinear method could possibly be more efficient in reducing the number of variables but would also have significant disadvantages. In particular, nonlinear methods typically require tuning specific parameters, objective criteria are often lacking, a proper weighting of observations is difficult, the methods are often not reversible, and it is harder to interpret the resulting indicators due to their nonlinear nature (Kraemer et al.2018). The salient feature of PCA is that an inverse projection is well defined and allows for a deeper inspection of the errors, which is not the case for nonlinear methods which learn a highly flexible transformation that is hard to invert. Therefore interpretability of the transform in meaningful physical units in the input space is often not possible. In the machine-learning community, this problem is known as the “pre-imaging problem” (Mika et al.1999; Arenas-Garcia et al.2013) and is a matter of current research.

4 Conclusions

To monitor the complexity of the changes occurring in times of an increasing human impact on the environment, we used PCA to construct indicators from a large number of data streams that track ecosystem state in space and time on a global scale. We showed that a large part of the variability of the terrestrial biosphere can be summarized using three indicators. The first emerging indicator represents carbon exchange, the second indicator shows the availability of water in the ecosystem, while the third indicator mostly represents a binary variable that indicates the presence of snow cover. The distribution in the space of the first two principal components reflects the general limitations of ecosystem productivity. Ecosystem production can be limited by either water or energy.

The first three indicators can detect many well-known phenomena without analyzing variables separately due to their compound nature. We showed that the indicators are capable of detecting seasonal hysteresis effects in ecosystems, as well as breakpoints, e.g. large-scale deforestation. The indicators can also track other changes to the seasonal cycle such as patterns of changes to the seasonal amplitudes and trends in ecosystems. Deviations from the mean seasonal cycle of the trajectories indicate extreme events such as the large-scale droughts in the Amazon during 2005 and 2010 and the Russian heat wave of 2010. The events are detected in a similar fashion as with classical multivariate anomaly detection methods while directly providing information on the underlying variables.

Using multivariate indicators, we gain a high level overview of phenomena in ecosystems, and the method therefore provides an interesting tool for analyses where it is required to capture a wide range of phenomena which are not necessarily known a priori. Future research should consider nonlinearities, adding data streams describing other important biosphere variables (e.g. related to biodiversity and habitat quality), and including different subsystems, such as the atmosphere or the anthroposphere.

Appendix A: Description of variables

Variables used describing the biosphere can be found in Table 1. Here we provide a more complete description of all variables.

Black-sky albedo is the reflected fraction of total incoming radiation under direct hemispherical reflectance, i.e., direct illumination (Muller et al.2011). This dataset is the broadband surface albedo including the visible, the near-infrared, and the shortwave-infrared spectrum (400–3000 nm). It is derived from the SPOT4-VEGETATION, SPOT5-VEGETATION2, and MERIS satellite sensors.

White-sky albedo is the reflected fraction of total incoming radiation under bihemispherical reflectance, i.e., diffuse illumination (Muller et al.2011). Together with black-sky albedo it can be used to estimate the albedo under different illumination conditions. This dataset is the broadband surface albedo including the visible, the near-infrared, and the shortwave-infrared spectrum (400–3000 nm). This dataset is derived from the SPOT4-VEGETATION, SPOT5-VEGETATION2, and MERIS satellite sensors.

Evaporation (mm d−1) is the amount of water evaporated per day, depending on the amount of available water and energy. This dataset is based on the GLEAMv3 model (Martens et al.2017), using satellite data from ESA CCI and SMOS to derive a number of variables.

Evaporative stress is modeled water stress for plants. Zero means that the vegetation has no water available for transpiration and 1 means that transpiration equals potential transpiration. This dataset is based on the GLEAMv3 model (Martens et al.2017), using satellite data from ESA CCI and SMOS to derive a number of variables.

fAPAR is the fraction of absorbed photosynthetically active radiation, a proxy for plant productivity (Disney et al.2016). This dataset is based on the GlobAlbedo dataset (, last access: 23 April 2020) and the MODIS fAPAR and leaf area index (LAI) products.

Gross primary productivity (GPP) is (gC m−2 d−1) the total amount of carbon fixed by photosynthesis (Tramontana et al.2016). This dataset is derived from upscaling eddy covariance tower observations to a global scale using machine-learning methods.

Terrestrial ecosystem respiration (TER) is (gC m−2 d−1) the total amount of carbon respired by the ecosystem, including autotrophic and heterotrophic respiration (Tramontana et al.2016). This dataset is derived from upscaling eddy covariance tower observations to a global scale using machine-learning methods.

Net ecosystem exchange (NEE) is (gC m−2 d−1) the total exchange of carbon of the ecosystem with the atmosphere NEE=GPP-TER (Tramontana et al.2016). This dataset is derived from upscaling eddy covariance tower observations to a global scale using machine-learning methods.

Latent energy (LE) is (W m−2) the amount of energy lost by the surface due to evaporation (Tramontana et al.2016). This dataset is derived from upscaling eddy covariance tower observations to a global scale using machine-learning methods.

Sensible heat (H) is (W m−2) the amount of energy lost by the surface due to radiation (Tramontana et al.2016). This dataset is derived from upscaling eddy covariance tower observations to a global scale using machine-learning methods.

Root-zone soil moisture is (m3 m−3) the moisture content of the root zone. This dataset is based on the GLEAMv3 model (Martens et al.2017), using satellite data from ESA CCI and SMOS to derive a number of variables.

Surface soil moisture is (mm3 mm−3) the soil moisture content at the soil surface. This dataset is based on the GLEAMv3 model (Martens et al.2017), using satellite data from ESA CCI and SMOS to derive a number of variables.

Appendix B: Time–space patterns of Components 1–3

Figure B1Time and space patterns of PC1–PC3, where the cut points are the same as in Fig. 8. The brown–green contrast shows the state of PC1, from low to high productivity. The blue–red contrast shows the state of PC2, from cold to dry. The brown–purple contrast shows the state of PC3, from dark to light. Panels (a), (e), and (i) are maps showing the state of PC1–PC3, respectively, on 1 January 2001. Panels (b), (c), and (d) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (a). Panels (f), (g), and (h) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (e). Panels (j), (k), and (l) show longitudinal cuts of PC1–PC3, respectively, at the red vertical line in (i).

Appendix C: Mean seasonal cycle extrema

Figure C1The minimum (a, c, e) and maximum (b, d, f) mean seasonal cycles of GPP (a, b), latent heat (c, d), and sensible heat (e, f). This illustrates the similarity of possibly very different ecosystems in terms of productivity and limitations. During peak growing season, many midlatitude areas have a similar productivity and latent energy release as tropical rain forests (b, d). The highest maximum seasonal sensible heat loss can be found in dry areas around the world and is lowest in areas with a wet climate such as tropical rain forests and maritime climates (f).

Appendix D: Spatial covariances of the components

Figure D1Pairwise covariances of the first three principal components mean seasonal cycles by space. (a) cov(PC1,PC2), (b) cov(PC1,PC3), and (c) cov(PC2,PC3). The bar charts show the distribution of the covariances. It can be seen that although two principal components are globally uncorrelated by their way of construction, they covary locally.

Appendix E: Changes in the seasonal amplitude

Figure E1Trends in the amplitude of the yearly cycle, 2001–2011. Only Theil–Sen estimators for significant slopes (p<0.05, unadjusted) are shown. Because there is only a single amplitude per year and therefore only 11 data points per time series, the Benjamini–Hochberg adjusted p values are not significant.

Appendix F: Breakpoints in trajectories

Figure F1Breakpoint detection, (a) on PC1, (b) on PC2, and (c) on PC3. The color indicates the year of the biggest breakpoint if a significant breakpoint was found, with gray if there was no significant breakpoint found.

As the environmental conditions change, due to climate change and human intervention, the local ecosystems may change gradually or abruptly. Detecting these changes is very important for monitoring the impact of climate change and land use change on the ecosystems. We applied breakpoint detection to the trajectories (Fig. F1).

Breakpoints on the first component were found in the entire Amazon, and the largest breakpoint is dated to the year 2005 during the large drought event. The entire eastern part of Australia shows its largest breakpoint towards the end of the time series because of a La Niña event, which caused lower temperatures and higher rainfall than usual during the years 2010 and 2011.

Code and data availability

The data are available and can be processed at, last access: 30 March 2020. The exact dataset and a docker container to reproduce the analysis can be found under (Kraemer et al.2020). The code to reproduce this analysis is available under (Kraemer2020) and, last access: 23 April 2020.

Author contributions

GK and MDM designed the study in collaboration with MR and GCV. GK conducted the analysis and wrote the manuscript with contributions from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.


We thank Fabian Gans and German Poveda for useful discussions. We thank Jake Nelson for proofreading a previous version of the manuscript. We thank Gregory Duveiller and the three anonymous reviewers for very helpful suggestions and Kirsten Thonicke for editorial advice that improved the manuscript greatly.

Financial support

This study is funded by the Earth System Data Lab – a project by the European Space Agency. Miguel D. Mahecha and Markus Reichstein have been supported by the Horizon 2020 EU project BACI under grant agreement no. 640176. Gustau Camps-Valls' work has been supported by the EU under the ERC consolidator grant SEDAL-647423.

The article processing charges for this open-access
publication were covered by the Max Planck Society.

Review statement

This paper was edited by Kirsten Thonicke and reviewed by Gregory Duveiller and three anonymous referees.


Abatzoglou, J. T., Rupp, D. E., and Mote, P. W.: Seasonal Climate Variability and Change in the Pacific Northwest of the United States, J. Clim., 27, 2125–2142,, 2014. a

Anav, A., Friedlingstein, P., Beer, C., Ciais, P., Harper, A., Jones, C., Murray-Tortarolo, G., Papale, D., Parazoo, N. C., Peylin, P., Piao, S., Sitch, S., Viovy, N., Wiltshire, A., and Zhao, M.: Spatiotemporal patterns of terrestrial gross primary production: A review: GPP Spatiotemporal Patterns, Rev. Geophys., 53, 785–818,, 2015. a, b

Aragão, L. E. O. C., Anderson, L. O., Fonseca, M. G., Rosan, T. M., Vedovato, L. B., Wagner, F. H., Silva, C. V. J., Silva Junior, C. H. L., Arai, E., Aguiar, A. P., Barlow, J., Berenguer, E., Deeter, M. N., Domingues, L. G., Gatti, L., Gloor, M., Malhi, Y., Marengo, J. A., Miller, J. B., Phillips, O. L., and Saatchi, S.: 21st Century Drought-Related Fires Counteract the Decline of Amazon Deforestation Carbon Emissions, Nat. Commun., 9, 146–149,, 2018. a

Ardisson, P.-L., Bourget, E., and Legendre, P.: Multivariate Approach to Study Species Assemblages at Large Spatiotemporal Scales: The Community Structure of the Epibenthic Fauna of the Estuary and Gulf of St. Lawrence, Can. J. Fish. Aquat. Sci., 47, 1364–1377,, 1990. a

Arenas-Garcia, J., Petersen, K. B., Camps-Valls, G., and Hansen, L. K.: Kernel Multivariate Analysis Framework for Supervised Subspace Learning: A Tutorial on Linear and Kernel Multivariate Methods, IEEE Signal Processing Magazine, 30, 16–29,, 2013. a

Babst, F., Poulter, B., Bodesheim, P., Mahecha, M. D., and Frank, D. C.: Improved tree-ring archives will support earth-system science, Nat. Ecol. Evol., 1, 1–2, 2017. a

Baldocchi, D. D.: How Eddy Covariance Flux Measurements Have Contributed to Our Understanding of Global Change Biology, Glob. Change Biol., 26, 242–260,, 2020. a, b

Barriopedro, D., Fischer, E. M., Luterbacher, J., Trigo, R. M., and García-Herrera, R.: The Hot Summer of 2010: Redrawing the Temperature Record Map of Europe, Science, 332, 220–224,, 2011. a

Beisner, B., Haydon, D., and Cuddington, K.: Alternative Stable States in Ecology, Front. Ecol. Environ., 1, 376–382,[0376:ASSIE]2.0.CO;2, 2003. a

Benjamini, Y. and Hochberg, Y.: Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing, J. Roy. Stat. Soc. B, 57, 289–300, 1995. a

Berger, M., Moreno, J., Johannessen, J. A., Levelt, P. F., and Hanssen, R. F.: ESA's Sentinel Missions in Support of Earth System Science, Remote Sens. Environ., 120, 84–90,, 2012. a

Blonder, B., Moulton, D. E., Blois, J., Enquist, B. J., Graae, B. J., Macias-Fauria, M., McGill, B., Nogué, S., Ordonez, A., Sandel, B., and Svenning, J.-C.: Predictability in Community Dynamics, Ecol. Lett., 20, 293–306,, 2017. a

Bowen, I. S.: The Ratio of Heat Losses by Conduction and by Evaporation from Any Water Surface, Phys. Rev., 27, 779–787,, 1926. a, b

Brown, R. L., Durbin, J., and Evans, J. M.: Techniques for Testing the.ournal of the Roy. Stat. Soc. B, 37, 149–192, 1975. a

Cattell, R. B.: The Scree Test For The Number Of Factors, Multivar. Behav. Res., 1, 245–276,, 1966. a

Chapin, F. S., Woodwell, G. M., Randerson, J. T., Rastetter, E. B., Lovett, G. M., Baldocchi, D. D., Clark, D. A., Harmon, M. E., Schimel, D. S., Valentini, R., Wirth, C., Aber, J. D., Cole, J. J., Goulden, M. L., Harden, J. W., Heimann, M., Howarth, R. W., Matson, P. A., McGuire, A. D., Melillo, J. M., Mooney, H. A., Neff, J. C., Houghton, R. A., Pace, M. L., Ryan, M. G., Running, S. W., Sala, O. E., Schlesinger, W. H., and Schulze, E.-D.: Reconciling Carbon-Cycle Concepts, Terminol. Method. Ecosys., 9, 1041–1050,, 2006. a

Chen, C., Park, T., Wang, X., Piao, S., Xu, B., Chaturvedi, R. K., Fuchs, R., Brovkin, V., Ciais, P., Fensholt, R., Tømmervik, H., Bala, G., Zhu, Z., Nemani, R. R., and Myneni, R. B.: China and India Lead in Greening of the World through Land-Use Management, Nature Sustainability, 2, 122–129,, 2019. a

Ciais, P., Reichstein, M., Viovy, N., Granier, A., Ogée, J., Allard, V., Aubinet, M., Buchmann, N., Bernhofer, C., Carrara, A., Chevallier, F., Noblet, N. D., Friend, A. D., Friedlingstein, P., Grünwald, T., Heinesch, B., Keronen, P., Knohl, A., Krinner, G., Loustau, D., Manca, G., Matteucci, G., Miglietta, F., Ourcival, J. M., Papale, D., Pilegaard, K., Rambal, S., Seufert, G., Soussana, J. F., Sanz, M. J., Schulze, E. D., Vesala, T., and Valentini, R.: Europe-Wide Reduction in Primary Productivity Caused by the Heat and Drought in 2003, Nature, 437, 529–533,, 2005. a

de Jong, R., de Bruin, S., de Wit, A., Schaepman, M. E., and Dent, D. L.: Analysis of Monotonic Greening and Browning Trends from Global NDVI Time-Series, Remote Sens. Environ., 115, 692–702,, 2011. a, b

Díaz, S., Settele, J., Brondízio, E. S., Ngo, H. T., Agard, J., Arneth, A., Balvanera, P., Brauman, K. A., Butchart, S. H. M., Chan, K. M. A., Garibaldi, L. A., Ichii, K., Liu, J., Subramanian, S. M., Midgley, G. F., Miloslavich, P., Molnár, Z., Obura, D., Pfaff, A., Polasky, S., Purvis, A., Razzaque, J., Reyers, B., Chowdhury, R. R., Shin, Y.-J., Visseren-Hamakers, I., Willis, K. J., and Zayas, C. N.: Pervasive Human-Driven Decline of Life on Earth Points to the Need for Transformative Change, Science, 366, 6471,, 2019. a

Disney, M., Muller, J.-P., Kharbouche, S., Kaminski, T., Voßbeck, M., Lewis, P., and Pinty, B.: A New Global fAPAR and LAI Dataset Derived from Optimal Albedo Estimates: Comparison with MODIS Products, Remote Sens., 8, 1–29,, 2016. a, b

Doughty, C. E., Metcalfe, D. B., Girardin, C. a. J., Amézquita, F. F., Cabrera, D. G., Huasco, W. H., Silva-Espejo, J. E., Araujo-Murakami, A., da Costa, M. C., Rocha, W., Feldpausch, T. R., Mendoza, A. L. M., da Costa, A. C. L., Meir, P., Phillips, O. L., and Malhi, Y.: Drought impact on forest carbon dynamics and fluxes in Amazonia, Nature, 519, 78–82,, 2015. a

Feldpausch, T. R., Phillips, O. L., Brienen, R. J. W., Gloor, E., Lloyd, J., Lopez-Gonzalez, G., Monteagudo-Mendoza, A., Malhi, Y., Alarcón, A., Dávila, E. A., Alvarez-Loayza, P., Andrade, A., Aragao, L. E. O. C., Arroyo, L., C, G. A. A., Baker, T. R., Baraloto, C., Barroso, J., Bonal, D., Castro, W., Chama, V., Chave, J., Domingues, T. F., Fauset, S., Groot, N., Coronado, E. H., Laurance, S., Laurance, W. F., Lewis, S. L., Licona, J. C., Marimon, B. S., Marimon-Junior, B. H., Bautista, C. M., Neill, D. A., Oliveira, E. A., dos Santos, C. O., Camacho, N. C. P., Pardo-Molina, G., Prieto, A., Quesada, C. A., Ramírez, F., Ramírez-Angulo, H., Réjou-Méchain, M., Rudas, A., Saiz, G., Salomão, R. P., Silva-Espejo, J. E., Silveira, M., ter Steege, H., Stropp, J., Terborgh, J., Thomas-Caesar, R., van der Heijden, G. M. F., Martinez, R. V., Vilanova, E., and Vos, V. A.: Amazon Forest Response to Repeated Droughts, Global Biogeochem. Cy., 30, 964–982,, 2016. a

Flach, M., Gans, F., Brenning, A., Denzler, J., Reichstein, M., Rodner, E., Bathiany, S., Bodesheim, P., Guanche, Y., Sippel, S., and Mahecha, M. D.: Multivariate anomaly detection for Earth observations: a comparison of algorithms and feature extraction techniques, Earth Syst. Dynam., 8, 677–696,, 2017. a

Flach, M., Sippel, S., Gans, F., Bastos, A., Brenning, A., Reichstein, M., and Mahecha, M. D.: Contrasting biosphere responses to hydrometeorological extremes: revisiting the 2010 western Russian heatwave, Biogeosciences, 15, 6067–6085,, 2018. a, b, c, d

Folke, C., Carpenter, S., Walker, B., Scheffer, M., Elmqvist, T., Gunderson, L., and Holling, C.: Regime Shifts, Resilience, and Biodiversity in Ecosystem Management, Annu. Rev. Ecol. Evol. S., 35, 557–581,, 2004. a

Forkel, M., Carvalhais, N., Verbesselt, J., Mahecha, M., Neigh, C., Reichstein, M., Forkel, M., Carvalhais, N., Verbesselt, J., Mahecha, M. D., Neigh, C. S. R., and Reichstein, M.: Trend Change Detection in NDVI Time Series: Effects of Inter-Annual Variability and Methodology, Remote Sens., 5, 2113–2144,, 2013. a

Forkel, M., Migliavacca, M., Thonicke, K., Reichstein, M., Schaphoff, S., Weber, U., and Carvalhais, N.: Codominant Water Control on Global Interannual Variability and Trends in Land Surface Phenology and Greenness, Glob. Change Biol., 21, 3414–3435,, 2015. a

Forkel, M., Carvalhais, N., Rodenbeck, C., Keeling, R., Heimann, M., Thonicke, K., Zaehle, S., and Reichstein, M.: Enhanced Seasonal CO2 Exchange Caused by Amplified Plant Productivity in Northern Ecosystems, Science, 351, 696–699,, 2016. a, b

Foster, A. C., Armstrong, A. H., Shuman, J. K., Shugart, H. H., Rogers, B. M., Mack, M. C., Goetz, S. J., and Ranson, K. J.: Importance of Tree- and Species-Level Interactions with Wildfire, Climate, and Soils in Interior Alaska: Implications for Forest Change under a Warming Climate, Ecol. Modell., 409, 108765,, 2019. a

García-Herrera, R., Díaz, J., Trigo, R. M., Luterbacher, J., and Fischer, E. M.: A Review of the European Summer Heat Wave of 2003, Crit. Rev. Env. Sci. Tec., 40, 267–306,, 2010. a

Gaumont-Guay, D., Black, T. A., Griffis, T. J., Barr, A. G., Jassal, R. S., and Nesic, Z.: Interpreting the Dependence of Soil Respiration on Soil Temperature and Water Content in a Boreal Aspen Stand, Agr. Forest Meteorol., 140, 220–235,, 2006. a

Graven, H. D., Keeling, R. F., Piper, S. C., Patra, P. K., Stephens, B. B., Wofsy, S. C., Welp, L. R., Sweeney, C., Tans, P. P., Kelley, J. J., Daube, B. C., Kort, E. A., Santoni, G. W., and Bent, J. D.: Enhanced Seasonal Exchange of CO2 by Northern Ecosystems Since 1960, Science, 341, 1085–1089,, 2013. a

Hao, Z., Zheng, J., Ge, Q., and Wang, W.: Historical Analogues of the 2008 Extreme Snow Event over Central and Southern China, Clim. Res., 50, 161–170,, 2011. a

Hendon, H. H., Lim, E.-P., Arblaster, J. M., and Anderson, D. L. T.: Causes and Predictability of the Record Wet East Australian Spring 2010, Clim. Dynam., 42, 1155–1174,, 2014. a

Higham, N. J.: The Accuracy of Floating Point Summation, SIAM J. Sci. Comput., 14, 783–799,, 1993. a

Horridge, M., Madden, J., and Wittwer, G.: The Impact of the 2002–2003 Drought on Australia, J. Policy Model., 27, 285–308,, 2005. a

Huang, K., Xia, J., Wang, Y., Ahlström, A., Chen, J., Cook, R. B., Cui, E., Fang, Y., Fisher, J. B., Huntzinger, D. N., Li, Z., Michalak, A. M., Qiao, Y., Schaefer, K., Schwalm, C., Wang, J., Wei, Y., Xu, X., Yan, L., Bian, C., and Luo, Y.: Enhanced Peak Growth of Global Vegetation and Its Key Mechanisms, Nat. Ecol. Evol., 2, 1897–1905,, 2018. a

IPBES: Summary for Policymakers of the Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, Summary for policymakers, IPBES, 39 pp., 2019. a

IPCC: Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Tech. rep., IPCC, Geneva, Swizerland, 2014. a

Ivits, E., Horion, S., Fensholt, R., and Cherlet, M.: Drought Footprint on European Ecosystems between 1999 and 2010 Assessed by Remotely Sensed Vegetation Phenology and Productivity, Glob. Change Biol., 20, 581–593,, 2014. a

Jolly, W. M., Cochrane, M. A., Freeborn, P. H., Holden, Z. A., Brown, T. J., Williamson, G. J., and Bowman, D. M. J. S.: Climate-Induced Variations in Global Wildfire Danger from 1979 to 2013, Nat. Commun., 6, 7537,, 2015. a

Jung, M., Koirala, S., Weber, U., Ichii, K., Gans, F., Camps-Valls, G., Papale, D., Schwalm, C., Tramontana, G., and Reichstein, M.: The FLUXCOM Ensemble of Global Land-Atmosphere Energy Fluxes, Sci. Data, 6, 1–14,, 2019. a, b, c, d, e, f, g

Keeling, C. D., Chin, J. F. S., and Whorf, T. P.: Increased Activity of Northern Vegetation Inferred from AtmosphericCO2 Measurements, Nature, 382, 146–149,, 1996. a

Kendall, M. G.: Rank Correlation Methods, Griffin, London, 202 pp., 1970. a

Khanna, J., Medvigy, D., Fueglistaler, S., and Walko, R.: Regional dry-season climate changes due to three decades of Amazonian deforestation, Nat. Clim. Change, 7, 200–204,, 2017. a

Kottek, M., Grieser, J., Beck, C., Rudolf, B., and Rubel, F.: World Map of the Köppen-Geiger climate classification updated, Meteorol. Z., 15, 259–263,, a

Kraemer, G.: gdkrmr/summarizing_the_state_of_the_biosphere v1.1.1, Zenodo,, 2020. a

Kraemer, G., Reichstein, M., and Mahecha, M. D.: dimRed and coRanking – Unifying Dimensionality Reduction in R, R J., 10, 342–358,, 2018. a, b

Kraemer, G., Camps-Valls, G., Reichstein, M., and Mahecha, M. D.: Summarizing the state of the terrestrial biosphere in few dimensions, Zenodo,, 2020. a

Kuan, C.-M. and Hornik, K.: The Generalized Fluctuation Test: A Unifying View, Economet. Rev., 14, 135–161,, 1995. a

Legendre, P. and Legendre, L.: Numerical Ecology: Second English Edition, Dev. Environ. Model., 20, 852 pp., 1998. a, b

Legendre, P., Planas, D., and Auclair, M.-J.: Succession des communautés de gastéropodes dans deux milieux différant par leur degré d'eutrophisation, Can. J. Zool., 62, 2317–2327,, 1984. a

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping Elements in the Earth's Climate System, P. Natl. Acad. Sci. USA, 105, 1786–1793,, 2008. a

Mahecha, M. D., Martínez, A., Lischeid, G., and Beck, E.: Nonlinear Dimensionality Reduction: Alternative Ordination Approaches for Extracting and Visualizing Biodiversity Patterns in Tropical Montane Forest Vegetation Data, Ecol. Inform., 2, 138–149,, 2007a. a, b

Mahecha, M. D., Reichstein, M., Lange, H., Carvalhais, N., Bernhofer, C., Grünwald, T., Papale, D., and Seufert, G.: Characterizing Ecosystem-Atmosphere Interactions from Short to Interannual Time Scales, Biogeosciences, 4, 743–758,, 2007b. a

Mahecha, M. D., Gans, F., Sippel, S., Donges, J. F., Kaminski, T., Metzger, S., Migliavacca, M., Papale, D., Rammig, A., and Zscheischler, J.: Detecting Impacts of Extreme Events with Ecological in Situ Monitoring Networks, Biogeosciences, 14, 4255–4277,, 2017. a

Mahecha, M. D., Gans, F., Brandt, G., Christiansen, R., Cornell, S. E., Fomferra, N., Kraemer, G., Peters, J., Bodesheim, P., Camps-Valls, G., Donges, J. F., Dorigo, W., Estupinan-Suarez, L. M., Gutierrez-Velez, V. H., Gutwin, M., Jung, M., Londoño, M. C., Miralles, D. G., Papastefanou, P., and Reichstein, M.: Earth system data cubes unravel global multivariate dynamics, Earth Syst. Dynam., 11, 201–234,, 2020. a

Mann, H. B.: Nonparametric Tests Against Trend, Econometrica, 13, 245–259,, 1945. a

Martens, B., Miralles, D. G., Lievens, H., Schalie, R. v. d., Jeu, R. A. M. d., Fernández-Prieto, D., Beck, H. E., Dorigo, W. A., and Verhoest, N. E. C.: GLEAM v3: satellite-based land evaporation and root-zone soil moisture, Geosci. Model Dev., 10, 1903–1925,, 2017. a, b, c, d, e, f, g, h, i

Metzger, M. J., Bunce, R. G. H., Jongman, R. H. G., Sayre, R., Trabucco, A., and Zomer, R.: A High-resolution Bioclimate Map of the World: A Unifying Framework for Global Biodiversity Research and Monitoring, Glob. Ecol. Biogeogr., 22, 630–638,, 2013. a

Mika, S., Scholkopf, B., Smola, A., Muller, K., Scholz, M., and Ratsch, G.: Kernel PCA and De-Noising in Feature Spaces, in: Advances In Neural Information Processing Systems, edited by: Kearns, M. S., Solla, S. A., and Cohn, D. A., Vol. 11 of Advances in Neural Information Processing Systems, 12th Annual Conference on Neural Information Processing Systems (NIPS), Denver, CO, 30 November–5 December 1998, 536–542, 1999. a

Miralles, D. G., Teuling, A. J., van Heerwaarden, C. C., and Vilà-Guerau de Arellano, J.: Mega-heatwave temperatures due to combined soil desiccation and atmospheric heat accumulation, Nat. Geosci., 7, 345–349,, 2014. a

Muller, J.-P., Lewis, P., Fischer, J., North, P., and Framer, U.: The ESA GlobAlbedo Project for mapping the Earth’s land surface albedo for 15 years from European sensors, Geophys. Res. Abstr., 13, EGU2011-10969, 2011. a, b, c, d

Najafi, E., Pal, I., and Khanbilvardi, R.: Climate Drives Variability and Joint Variability of Global Crop Yields, Sci. Total Environ.t, 662, 361–372,, 2019. a

Nasahara, K. N. and Nagai, S.: Review: Development of an in Situ Observation Network for Terrestrial Ecological Remote Sensing: The Phenological Eyes Network (PEN), Ecol. Res., 30, 211–223,, 2015. a

Nicholls, N.: The Changing Nature of Australian Droughts, Climatic Change, 63, 323–336,, 2004. a

Nicholson, S. E.: A detailed look at the recent drought situation in the Greater Horn of Africa, J. Arid Environ., 103, 71–79,, 2014. a

Papale, D., Black, T. A., Carvalhais, N., Cescatti, A., Chen, J., Jung, M., Kiely, G., Lasslop, G., Mahecha, M. D., Margolis, H., Merbold, L., Montagnani, L., Moors, E., Olesen, J. E., Reichstein, M., Tramontana, G., van Gorsel, E., Wohlfahrt, G., and Ráduly, B.: Effect of Spatial Sampling from European Flux Towers for Estimating Carbon and Water Fluxes with Artificial Neural Networks, J. Geophys. Res.-Biogeo., 120, 1941–1957,, 2015. a

Parmesan, C.: Ecological and Evolutionary Responses to Recent Climate Change, Ann. Rev. Ecol. Evol. S., 37, 637–669,, 2006. a

Pearson, K.: On Lines and Planes of Closest Fit to Systems of Points in Space, Philos. Mag., 2, 559–572, 1901. a

Piao, S., Wang, X., Park, T., Chen, C., Lian, X., He, Y., Bjerke, J. W., Chen, A., Ciais, P., Tømmervik, H., Nemani, R. R., and Myneni, R. B.: Characteristics, drivers and feedbacks of global greening, Nat. Rev. Earth Environ., 1, 14–27,, 2019. a

Piao, S., Wang, X., Wang, K., Li, X., Bastos, A., Canadell, J. G., Ciais, P., Friedlingstein, P., and Sitch, S.: Interannual Variation of Terrestrial Carbon Cycle: Issues and Perspectives, Glob. Change Biol., 26, 300–318,, 2020. a

Rao, M., Saw Htun, Platt, S. G., Tizard, R., Poole, C., Than Myint, and Watson, J. E. M.: Biodiversity Conservation in a Changing Climate: A Review of Threats and Implications for Conservation Planning in Myanmar, AMBIO, 42, 789–804,, 2013. a

Reichstein, M., Bahn, M., Ciais, P., Frank, D., Mahecha, M. D., Seneviratne, S. I., Zscheischler, J., Beer, C., Buchmann, N., Frank, D. C., Papale, D., Rammig, A., Smith, P., Thonicke, K., van der Velde, M., Vicca, S., Walz, A., and Wattenbach, M.: Climate extremes and the carbon cycle, Nature, 500, 287–295,, 2013. a

Renner, M., Brenner, C., Mallick, K., Wizemann, H.-D., Conte, L., Trebs, I., Wei, J., Wulfmeyer, V., Schulz, K., and Kleidon, A.: Using Phase Lags to Evaluate Model Biases in Simulating the Diurnal Cycle of Evapotranspiration: A Case Study in Luxembourg, Hydrol. Earth Syst. Sci., 23, 515–535,, 2019. a

Richardson, A. D., Braswell, B. H., Hollinger, D. Y., Burman, P., Davidson, E. A., Evans, R. S., Flanagan, L. B., Munger, J. W., Savage, K., Urbanski, S. P., and Wofsy, S. C.: Comparing Simple Respiration Models for Eddy Flux and Dynamic Chamber Data, Agr. Forest Meteorol., 141, 219–234,, 2006. a

Rosenfeld, D., Zhu, Y., Wang, M., Zheng, Y., Goren, T., and Yu, S.: Aerosol-Driven Droplet Concentrations Dominate Coverage and Water of Oceanic Low-Level Clouds, Science, 363, eaav0566,, 2019. a

Sarmah, S., Jia, G., and Zhang, A.: Satellite View of Seasonal Greenness Trends and Controls in South Asia, Environ. Res. Lett., 13, 034026,, 2018. a

Schimel, D. and Schneider, F. D.: Flux Towers in the Sky: Global Ecology from Space, New Phytol., 224, 570–584,, 2019. a

Schwartz, M. D.: Monitoring Global Change with Phenology: The Case of the Spring Green Wave, Int. J. Biometeorol., 38, 18–22,, 1994. a

Schwartz, M. D.: Green-Wave Phenology, Nature, 394, 839–840,, 1998. a, b

Sen, P. K.: Estimates of the Regression Coefficient Based on Kendall's Tau, J. Am. Stat. Assoc., 63, 1379–1389,, 1968. a

Sippel, S., Reichstein, M., Ma, X., Mahecha, M. D., Lange, H., Flach, M., and Frank, D.: Drought, Heat, and the Carbon Cycle: A Review, Current Climate Change Reports, 4, 266–286,, 2018. a

Sitch, S., Friedlingstein, P., Gruber, N., Jones, S. D., Murray-Tortarolo, G., Ahlström, A., Doney, S. C., Graven, H., Heinze, C., Huntingford, C., Levis, S., Levy, P. E., Lomas, M., Poulter, B., Viovy, N., Zaehle, S., Zeng, N., Arneth, A., Bonan, G., Bopp, L., Canadell, J. G., Chevallier, F., Ciais, P., Ellis, R., Gloor, M., Peylin, P., Piao, S. L., Le Quéré, C., Smith, B., Zhu, Z., and Myneni, R.: Recent Trends and Drivers of Regional Sources and Sinks of Carbon Dioxide, Biogeosciences, 12, 653–679,, 2015. a

Song, X.-P., Hansen, M. C., Stehman, S. V., Potapov, P. V., Tyukavina, A., Vermote, E. F., and Townshend, J. R.: Global Land Change from 1982 to 2016, Nature, 560, 639–643,, 2018. a, b

Steffen, W., Richardson, K., Rockström, J., Cornell, S. E., Fetzer, I., Bennett, E. M., Biggs, R., Carpenter, S. R., de Vries, W., de Wit, C. A., Folke, C., Gerten, D., Heinke, J., Mace, G. M., Persson, L. M., Ramanathan, V., Reyers, B., and Sörlin, S.: Planetary Boundaries: Guiding Human Development on a Changing Planet, Science, 347, 1259855-1–1259855-10,, 2015. a

Stine, A. R., Huybers, P., and Fung, I. Y.: Changes in the Phase of the Annual Cycle of Surface Temperature, Nature, 457, 435–440,, 2009. a, b

Tang, J., Baldocchi, D. D., and Xu, L.: Tree Photosynthesis Modulates Soil Respiration on a Diurnal Time Scale, Glob. Change Biol., 11, 1298–1304,, 2005. a

Theil, H.: A Rank-Invariant Method of Linear and Polynomial Regression Analysis, I, II, III, Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, 53, 386–392, 1950a. a

Theil, H.: A Rank-Invariant Method of Linear and Polynomial Regression Analysis, I, II, III, Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, 53, 521–525, 1950b. a

Theil, H.: A Rank-Invariant Method of Linear and Polynomial Regression Analysis, I, II, III, Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, 53, 1397–1412, 1950c. a

Tramontana, G., Jung, M., Schwalm, C. R., Ichii, K., Camps-Valls, G., Ráduly, B., Reichstein, M., Arain, M. A., Cescatti, A., Kiely, G., Merbold, L., Serrano-Ortiz, P., Sickert, S., Wolf, S., and Papale, D.: Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms, Biogeosciences, 13, 4291–4313,, 2016. a, b, c, d, e, f, g, h, i, j, k

Van Der Maaten, L., Postma, E., and Van den Herik, J.: Dimensionality Reduction: A Comparative Review, J. Mach. Learn. Res., 10, 66–71, 2009. a

van der Maaten, L., Schmidtlein, S., and Mahecha, M. D.: Analyzing Floristic Inventories with Multiple Maps, Ecol. Inform., 9, 1–10,, 2012. a

Verbesselt, J., Hyndman, R., Newnham, G., and Culvenor, D.: Detecting Trend and Seasonal Changes in Satellite Image Time Series, Remote Sens. Environ., 114, 106–115,, 2010. a

Wilks, D. S.: Chapter 12 – Principal Component (EOF) Analysis, in: International Geophysics, edited by Wilks, D. S., vol. 100 of Statistical Methods in the Atmospheric Sciences, Academic Press, 519–562,, 2011. a

Wingate, L., Ogée, J., Cremonese, E., Filippa, G., Mizunuma, T., Migliavacca, M., Moisy, C., Wilkinson, M., Moureaux, C., Wohlfahrt, G., Hammerle, A., Hörtnagl, L., Gimeno, C., Porcar-Castell, A., Galvagno, M., Nakaji, T., Morison, J., Kolle, O., Knohl, A., Kutsch, W., Kolari, P., Nikinmaa, E., Ibrom, A., Gielen, B., Eugster, W., Balzarolo, M., Papale, D., Klumpp, K., Köstner, B., Grünwald, T., Joffre, R., Ourcival, J.-M., Hellstrom, M., Lindroth, A., George, C., Longdoz, B., Genty, B., Levula, J., Heinesch, B., Sprintsin, M., Yakir, D., Manise, T., Guyon, D., Ahrends, H., Plaza-Aguilar, A., Guan, J. H., and Grace, J.: Interpreting Canopy Development and Physiology Using a European Phenology Camera Network at Flux Sites, Biogeosciences, 12, 5995–6015,, 2015.  a

Wolter, K. and Timlin, M. S.: El Niño/Southern Oscillation Behaviour since 1871 as Diagnosed in an Extended Multivariate ENSO Index (MEI.Ext), Int. J. Climatol., 31, 1074–1087,, 2011. a

Yan, T., Song, H., Wang, Z., Teramoto, M., Wang, J., Liang, N., Ma, C., Sun, Z., Xi, Y., Li, L., and Peng, S.: Temperature Sensitivity of Soil Respiration across Multiple Time Scales in a Temperate Plantation Forest, Sci. Total Environ., 688, 479–485,, 2019. a

Zeileis, A., Leisch, F., Hornik, K., and Kleiber, C.: Strucchange: An R Package for Testing for Structural Change in Linear Regression Models, J. Stat. Softw., 7, 1–38,, 2002. a

Zeng, N., Zhao, F., Collatz, G. J., Kalnay, E., Salawitch, R. J., West, T. O., and Guanter, L.: Agricultural Green Revolution as a Driver of Increasing Atmospheric CO2 Seasonal Amplitude, Nature, 515, 394–397,, 2014. a

Zhang, Q., Phillips, R. P., Manzoni, S., Scott, R. L., Oishi, A. C., Finzi, A., Daly, E., Vargas, R., and Novick, K. A.: Changes in Photosynthesis and Soil Moisture Drive the Seasonal Soil Respiration-Temperature Hysteresis Relationship, Agr. Forest Meteorol., 259, 184–195,, 2018. a

Zhou, L., Tian, Y., Myneni, R. B., Ciais, P., Saatchi, S., Liu, Y. Y., Piao, S., Chen, H., Vermote, E. F., Song, C., and Hwang, T.: Widespread Decline of Congo Rainforest Greenness in the Past Decade, Nature, 509, 86–90,, 2014. a, b

Zhu, Z., Piao, S., Myneni, R. B., Huang, M., Zeng, Z., Canadell, J. G., Ciais, P., Sitch, S., Friedlingstein, P., Arneth, A., Cao, C., Cheng, L., Kato, E., Koven, C., Li, Y., Lian, X., Liu, Y., Liu, R., Mao, J., Pan, Y., Peng, S., Peñuelas, J., Poulter, B., Pugh, T. A. M., Stocker, B. D., Viovy, N., Wang, X., Wang, Y., Xiao, Z., Yang, H., Zaehle, S., and Zeng, N.: Greening of the Earth and Its Drivers, Nat. Clim. Change, 6, 791–795,, 2016. a, b, c, d, e

Publications Copernicus
Short summary
To closely monitor the state of our planet, we require systems that can monitor the observation of many different properties at the same time. We create indicators that resemble the behavior of many different simultaneous observations. We apply the method to create indicators representing the Earth's biosphere. The indicators show a productivity gradient and a water gradient. The resulting indicators can detect a large number of changes and extremes in the Earth system.
Final-revised paper