Introduction

Understanding how the future hydrological cycle will change at regional scales is crucial for impact assessment and adaptation planning1. Previous research has produced apparently conflicting insights on how land hydrological conditions have changed2,3,4,5, ranging from increasing aridity over land6 and further drying in the world’s driest regions7,8,9, to the acceleration of the hydrological cycle10 and increased vegetation “greening” (implying relaxation in water-limitations)11,12,13,14,15. Some studies have suggested that hydrological cycle trends are consistent with ‘wet gets wetter and dry gets drier’9,16, while others have found that both this and the opposite pattern ‘dry gets wetter and wet gets drier’ is detectable on land3,17. These conflicting findings may stem from different interpretations of ‘drying’, ‘aridity’ and ‘intensification’, depending on the metrics and spatial-scale aggregation methods used4,18. Changes in the terrestrial hydrological cycle have typically been explored by examining trends in just one3,5,9,10,19,20,21, or a subset16,22,23,24,25,26,27,28,29,30,31,32 of its components, or in aridity metrics that do not capture the full complexity of changes4. Furthermore, many studies use a single-source dataset, so that results are very sensitive to dataset choice2,4,24,33 and its inherent uncertainties4,21. Here, we undertake a simultaneous analysis of trends in all components of the hydrological cycle and vegetation greenness (based on Normalised Difference Vegetation Index; NDVI), using global observationally-derived gridded datasets of P, ET, Q and change in total water storage (ΔTWS) that are constrained with in situ measurements (Supplementary Table 1), over the period 1980–2012. Crucially, our approach imposes physical conservation constraints in calculating trends. The analysis of trends utilises multiple products and observational sources in the derivation of each budget component, as well as uncertainty bounds informed by observational constraints. Sampling within these uncertainty estimates was used to determine the robustness of the trend (see Methods). Changes in the hydrological cycle are described at the grid scale (0.25°) and over precipitation regimes derived by clustering the land based on precipitation mean and variability (Supplementary Fig. 1 and Supplementary Table 2).

Results

Global patterns of change in the hydrological cycle

We first examine the spatial distribution of robust trends in P, Q, ET, and NDVI over 1980–2012 (Fig. 1 – dark blue and red). We find that 51% of the land surface has experienced robust changes in at least one component of the hydrological cycle. P has changed robustly across ~13% of the land compared to ~34% and 20% for ET and Q, respectively. Regionally, P and Q show a broad tendency for increases in parts of the Sahel, southern Africa, Eurasia, the Amazon, southeast Asia and northern Australia. Declines in both P and Q are detected in southwestern North America, southern South America and parts of Asia. However, robust trends in P (dark colours in Fig. 1) are more limited in spatial extent and less spatially coherent than in Q. ET shows robust increases in most of Eurasia, southern Asia, and western North America. Declines in ET are detected in central Africa, western North America, the Middle East, and central South America. ET increases in northern mid- and high- latitudes and central South America are consistent with greening trends (i.e. an increase in maximum NDVI14, and a lengthening of the growing season). We also detect greening trends in some regions where no robust changes in the hydrological cycle are observed, particularly in Australia (Fig. 1).

Fig. 1: Spatial patterns of trends in P, Q, ET and NDVI.
figure 1

(Maps) Spatial patterns of trends in each of the annual totals in (a) P, (b) Q, and (c) ET, and (d) annual averages of monthly NDVI derived for the period 1980–2012 using Mann–Kendall and Sen’s slope methods. Grid cells in orange and light blue correspond to uncertain trends because the confidence interval of the slope encompasses a mix of negative and positive values. Grid cells in beige correspond to inconclusive trends either (i) because trend slopes computed for different estimates of the same component do not agree, (ii) due to artefacts in the data (NDVI), or (iii) because shorter period trends are inconsistent pre- and postdata assimilation technique. (Bars) Percentage of land area in each precipitation regime displaying these changes. The text on the top and bottom of each bar shows the areally-averaged slopes of positive and negative trends respectively. Precipitation regimes are labelled from driest to wettest as very dry with very high variability (V.dry, V.H.variability), dry with high variability (Dry, H.variability), dry with medium variability (Dry, M.variability), mild dry with medium variability (M.dry, M.variability), mild wet with medium variability (M.wet, M.variability), wet with medium variability (Wet, M.variability), wet with low variability (Wet, L.variability), very wet with low variability (V.wet, L.variability), and extremely wet with low variability (Ex.wet, L.variability). The spatial distribution of the precipitation regimes is illustrated in Supplementary Fig. 1.

The observed robust trends in P, ET and Q, while being less spatially widespread than previously reported changes in these hydroclimatic variables, do agree with many earlier studies on the overall spatial pattern of positive and negative trends2,32,34,35. We find notable disagreement, however, in Q trends with other reports22 particularly in central and southern South America and parts of Europe and India, likely stemming from different time periods in other work than in this study. Q trends in India closely resemble observed trends in regional studies36. Over most of the land, the directional changes in Q closely resemble those in P, despite our assessment being based on two independent datasets.

Trends in the change in the total water storage (ΔTWS) (Supplementary Fig. 2) are less spatially coherent than patterns in the other components of the hydrological cycle. Interestingly, these ΔTWS trends are not comparable to trends in total water storage (TWS)31, and instead indicate whether the rate of change in stored water equivalent is increasing or decreasing. ΔTWS shows robust increases in northern Australia, southern Sahel and some adjacent areas, and robust decreases in western China, in line with P trends. Robust declines are found in the western Amazon, parts of northern Siberia, and north-western Quebec, despite P increases. Elsewhere trends are highly localised and spatially variable.

Patterns of changes across precipitation regimes

We next examine trends aggregated by precipitation regimes to understand how trends vary from dry to wet environments. The magnitude of changes in P and Q consistently strengthen as we move from the very dry, very high variability to extremely wet, low variability precipitation regimes (Fig. 1- numbers above/below bar plots). The proportion of land exhibiting positive Q and P trends also increases from drier to wetter precipitation regimes and vice versa for negative trends. However, the proportion of land area showing increasing P and Q is much greater than the proportion of land area showing decreasing P and Q across most precipitation regimes, suggesting more widespread wetting than drying in all regimes. The trends are clearly spatially complex, and changes across multiple water balance components appear to be better explained by ‘wet gets wetter but dry does not get drier’ than ‘wet gets wetter and dry gets drier’9,16.

The magnitude of ET change (Fig. 1—numbers above/below bar plots) varies much less across precipitation regimes than P and Q, but shows the same tendency toward increasing magnitude from drier to wetter regimes. Unlike P and Q, decreasing ET is more common in the wet regimes than in the dry regimes and particularly in the wet, moderate variability regime. Similar to P and Q, increasing ET trends are more widespread than decreasing trends but robust increases cover a similar land area (~25%) across all precipitation regimes, with the exception of the wettest and driest class. The increasing trends in ET largely follow the pattern of NDVI increases but note NDVI increases are robust over a larger land area. Meanwhile, robust decreasing trends in ET are more widespread than for NDVI, particularly in the wetter regimes.

When grid cells are averaged across precipitation regimes, detected trends in P, ET, Q and NDVI are overwhelmingly positive (Fig. 2 and Supplementary Fig. 3); no statistically significant declines are detected in any regime or water balance component. These results remain true when trends in the dry precipitation regimes are computed separately at regions classified as non-drylands (‘Humid’ or ‘Cold’) and drylands (‘Hyper-arid’, ‘Arid’, ‘Semi-arid’, and ‘Dry Sub-humid’) based on common classification schemes (Supplementary Fig. 4). In the dry precipitation regimes, increases in the hydroclimatic variables mainly occur in non-drylands (Supplementary Fig. 3). This indicates that despite spatial variability within each precipitation regime, all the water balance components are increasing overall and the areas that exhibit positive trends dominate the trends. In each precipitation regime, positive trends in all hydroclimatic variables are either more widespread or of larger magnitudes than negative trends. These spatially aggregated results (Supplementary Fig. 2) therefore show that both ‘wet’ and ‘dry’ get wetter, although regional patterns remain complex (Fig. 1). The simultaneous increases in P, Q and ET identified across multiple precipitation regimes, and in particular the two wettest classes, also suggest a possible intensification of the water cycle as suggested in previous work10. ET, however, is changing at a much slower rate than P and Q likely due to the strong energy limitation of ET in these wet regions.

Fig. 2: Global trends in P, ET, Q and NDVI.
figure 2

Time series of annual totals in (a) P, (b) ET and (c) Q in (mm), and (d) annual averages of monthly NDVI (unitless × 10−4) computed for the dry and wet precipitation regimes during 1980–2012. Where trends are statistically significant, trendlines are displayed, and slopes and confidence intervals computed using Mann–Kendall and Sen’s slope methods are shown. Trends are deemed statistically significant if the confidence interval of the slope is strictly positive or strictly negative. The nine precipitation regimes are labelled from driest (bottom) to wettest (top) as listed in caption for Fig. 1.

Regional changes in the hydrological cycle and their likely drivers

Robust Q and ET trends are more widespread than those in P and we find many instances of unambiguous trends in ET and Q despite no robust trends in P (Fig. 1). Robust Q and ET trends occur across 44% of the regions where P trends are not robust (Fig. 3 middle row: sum of percentages represented in the (1) blue and (2) red ET bars, and (3) blue and red Q bars in the uncertain ET category), suggesting the trends in ET and Q are not merely driven by P changes as they can also occur due to changing atmospheric demand, vegetation responses to rising [CO2] or human intervention (i.e., land use and land cover change). Where positive P trends are observed (9.1% of the land), ~40% of these regions exhibit positive ET trends compared to ~6% of the land exhibiting negative ET trends (Fig. 3). This indicates that when P increases, ~87% of the robust trends in ET are also positive. We also find that where there is a negative P trend, ~73% of the robust trends in ET are also negative. Changes in Q are more strongly associated with P compared to ET. While only ~5% of the land area has robust trends in both Q and P, the proportion that shows a change in Q with the same directional change as P, is 95% when the change is positive, compared to 85% when it is negative.

Fig. 3: Spatial analysis of directional change in trends of P, ET, Q and NDVI.
figure 3

Percentage of land area displaying positive, negative, and uncertain trends in P (columns), and in each case, proportion of land displaying positive, negative, and uncertain trends in ET, and corresponding Q, and NDVI trends (rows). Uncertain trends include uncertain positive, uncertain negative and inconclusive trends.

While changes in ET typically follow changes in P, regional ET changes are largely associated with uncertain P trends (Fig. 4a, b) and in many regions, robust ET increases are associated with NDVI increases (e.g. Europe, India and China). These changes likely have a range of region-specific drivers, including greening due to increased plant water-use efficiency linked with rising [CO2]37, warming-induced changes in the growing season14, agricultural improvement in Europe15; and intensive farming and afforestation in India and southeast China15. In several regions, however, robust ET increases are associated with uncertain trends in both P and NDVI, including parts of Asia and North America (Fig. 4a) although significant P trends in eastern America have been identified when using a longer term dataset38. Past studies have suggested these changes may be related to increased vapour pressure deficit39 and increased radiation due to solar “brightening”40. The positive ET trend in northern high latitudes could also be associated with temperature-driven greening (i.e., a lengthening of the growing season)41,42 but poor NDVI data quality prevents the detection of robust trends in these regions. Robust decreases in ET are largely associated with declining or uncertain P and NDVI (Fig. 4b). Despite widespread ET trends, many challenges still exist in untangling the confounding effects of these drivers and other influences (e.g. temperature increases) on ET change, and such analysis is confounded by a lack of reliable, observational data in many regions.

Fig. 4: Spatial analysis of simultaneous change in trends of P, ET, Q and NDVI.
figure 4

Spatial distribution of positive (+), negative (−) and uncertain (u) trends in P and NDVI coincident with positive ET trends (a) and negative ET trends (b). Spatial distribution of positive, negative and uncertain trends in P and ET coincident with positive Q trends (c) and negative Q trends (d). Uncertain trends include uncertain positive, uncertain negative and disagreeing trends.

P and Q trends are both positive in parts of the Sahel, western Namibia, and northern Asia (Fig. 4c). In the Sahel, increasing precipitation has been attributed to anthropogenic climate change and is expected to continue in the future43. The observed changes in eastern Namibia are robust despite precipitation being highly variable in this region (very dry, very high variability regime) and are consistent across both P and Q. These results are inconsistent with some previous reports44 but agree with findings from a regional study45. In northern Asia, Q increases are likely caused by enhanced snow melting induced by climate change46. Positive P and Q trends are also observed in northeast Australia, and the Maritime Continent. Conversely, P and Q are both negative in the western US (Fig. 4d). These negative P and Q trends are coupled with decreasing ET which may relate to instances of droughts47. We detect negative Q trends in northern India despite no robust change in P, which is likely attributed to the regulation of water extraction for agriculture48. We also note that the grid-scale runoff trends are not robust over river basins in the Mediterranean and Northern Europe despite previous studies identifying significant basin-scale streamflow trends in these regions30,46. Spatial aggregation represents one potential source of this type of discrepancy in trend analyses.

Hydroclimatic trends may also highlight areas potentially more prone to flooding (Fig. 4c). For instance, both P and Q have increased in part of the Niger basin while ET has not changed, suggesting a possible increase in flood pre-conditioning, with other studies supporting the idea that this trend might continue49. P and Q have also increased in parts of Papua New Guinea and Indonesia, which are already prone to frequent and severe flooding. Similar trends in the Amazon agree with findings detecting recent intensification of flooding extremes in this region driven by a strengthening Walker circulation50.

We find that in the wet regions (e.g. Southeast Asia), changes in Q and P are closely related unlike changes in ET and P. In comparison, in the dry regions (e.g. in northern Africa), changes in ET and P are more alike than changes in Q and P. This finding confirms that in dry regions P is increasingly partitioned into ET rather than Q, and in wet regions P is increasingly partitioned into Q rather than ET51 with implications for freshwater availability in these water-limited regions.

Discussion

Robust changes in the hydrological cycle are observed over more than 50% of the global land area even when robust measures for assessing the significance of trends are implemented. Our trend analysis is enabled by an approach that uses observationally constrained components of the hydrological cycle from multiple sources (e.g. remote sensing, machine learning and land surface modelling) and considers their uncertainties and consistencies in a water balance conservation context. This conservative approach strengthens our confidence in the trends we observe. We find that ~38% of the land shows clear trends in ET or Q despite no robust changes in P. While ET is the most widely affected component of the hydrological cycle (~34% of the land), it is changing at a much slower rate than Q and P. Over large areas where trends in P are deemed uncertain, the confidence interval of the trends’ slopes includes a mix of positive and negative values. The inconclusive directional change in regional P is mostly attributed to strong connection of regional precipitation with the modes of climate variability. For example, the phase changes of the Pacific Decadal Oscillation (PDO) and El Niño-Southern Oscillation (ENSO) have been associated with marked changes in precipitation pattern anomalies worldwide, particularly in Australia, eastern Siberia and the Americas (Fig. 3 in Ref. 52). During the period of this analysis, a warm PDO phase occurred throughout 1980–1998 followed by a cool phase throughout 1999–2012. The transition between a warm to a cool PDO phase has likely contributed to the widespread inconclusive trends in P. On the other hand, there has been a slight tendency toward more La Niña than El Niño during this period, which might have contributed to the upward P trend in eastern Siberia and Northern South America, and downward P trend in Southwest North America and Central South America. Exploring changes in the hydrological cycle over longer time periods would likely help uncover new robust trends which have been inconclusive over 1980–2012. Global mean annual P, ET, and Q have increased at rates 0.8 mm × year−1, 0.3 mm × year−1, and 0.6 mm × year−1 respectively. These changes are equivalent to increases by 3%, 1.5% and 6% per 1°K of global mean surface warming for P, ET, and Q respectively.

Declining freshwater resources occur in regions that are already prone to water stress (e.g., the Middle East and western USA), in some intensive agricultural areas (e.g., key ‘breadbaskets’ in the Americas, Asia and Australia) and some densely populated areas including northern India and China. Negative trends in freshwater resources are expected to continue in most of these regions53, while water demand and/or population are expected to increase54 with implications for water scarcity in these regions55. These areas should be priority study regions for understanding the impact of changes in the hydrological cycle to global food security and to inform water policy regulation and risk, and adaptation management especially in low socioeconomic countries where vulnerability is high. We find that at least 6.3% of the world population is exposed to declining freshwater resources as evidenced by Q decreases (Supplementary Figure 5). In contrast, we detect a tendency toward increased freshwater resources in several other densely populated regions, including Indonesia, the Sahel and parts of southern Africa as evidenced by positive P and/or Q trends. Possible increase in flood pre-conditioning is identified in the Niger basin, Papua New Guinea, Indonesia, and the Amazon.

Our results demonstrate that hydroclimatic changes over land are multifaceted and not necessarily related to the background level of dryness or wetness of the land. In contrast to some earlier work9,16, our results support a general conclusion that over land “wet gets wetter but dry does not get drier”. Discrepancies between our work and previous studies may be explained by different choices in spatial aggregation, time period, variables used, employed datasets, and dryness definitions, making comparisons only partial. Here, for example, we use precipitation regimes, whereas other studies have used different spatial aggregation approaches, based on dryness (e.g. aridity indices) and/or spatial (e.g. river basin) characteristics51. Clustering land by precipitation regimes helps examine how the relationships between trends in P, ET, and Q play out in different environments, and there are numerous ways we could have chosen to separate these. Clustering simply as a function of the total annual precipitation and its seasonal variability allowed us to provide an intuitive separation, while minimising the number of clusters, so that we could explicitly explore the nature of trends and changes within each cluster, thereby covering the entire land surface. An alternative choice may have been to separate the dry precipitation regions into energy and water-limited clusters, which would result in a more natural latitudinal separation. We did not find that such a separation changed our conclusions (see Supplementary Fig. 3) and therefore opted not to pursue such an approach.

Understanding whether these observed trends will continue in the future remains a pressing challenge. The approach we use here to constrain hydrological cycle estimates is easily updated as new data sources become available. Next-generation of satellite missions (e.g. BIOMASS) may offer new understanding on the drivers and the scale of change (e.g. ECOSTRESS). Most importantly, by examining all components of the hydrological cycle, enforcing conservation, utilising observationally constrained uncertainty estimates and examining the statistical robustness of trends in this context, we now have an approach that can reconcile the apparent discrepancies in reported historical changes in the hydrological cycle, enabling new understanding that will help constrain future uncertainty.

Methods

Employed datasets

Hybrid datasets synthesised from multiple global estimates and constrained with in-situ measurement are used to analyse changes in the hydrological cycle over the global land for the period 1980–2012. These are MSWEP (V2.4)34, DOLCE (V3)56, and LORA (V1)57 for precipitation (P), runoff (Q), and evapotranspiration (ET) respectively. DOLCE is derived by merging four global ET datasets including physical-, empirical- and reanalysis-based estimates, and in-situ constraining ET data from more than 250 flux sites. LORA is derived by merging 11 global runoff datasets and in situ constraining streamflow records. Both DOLCE and LORA provide temporally and spatially variant uncertainty estimates. In DOLCE, uncertainty reflects the actual deviation from the measured ET at site locations58, while in LORA, uncertainty reflects how well aggregated runoff fields within a river basin reproduce streamflow in that basin59. MSWEP is derived by merging nine gauge-, satellite-, and reanalysis-based gridded datasets, and constraining in situ data from more than 7500 gauges34. The change in the total water storage (ΔTWS) is then taken as the residual of the surface water budget and computed by subtracting ET and Q from P, and used to analyse changes in ΔTWS during 1980–2012.

NDVI data from GIMMS (V1.1)60 is used to analyse changes in vegetation greening for the period 1982–2012. We resample the original bi‐monthly NDVI data to monthly time steps by taking the maximum of the bi‐monthly values for each calendar month to remove artefacts e.g., due to cloud contamination. We then mask pixels with < 85% good quality monthly values (i.e. quality flag = 0; Ref. 61).

Two additional datasets are incorporated to constrain the analysis of trends. These are the rain gauge-based precipitation dataset, GPCC62, and monthly ΔTWS derived from the monthly equivalent water thickness of water storage anomalies provided in GRACE Mascons dataset63 (denoted hereafter as ΔTWS-GRACE). The derivation of ΔTWS-GRACE follows the methods described in Ref. 64. As explained in ‘Section Data assimilation’, ΔTWS-GRACE is used as an input in a data assimilation technique to detect inconclusive trends in each component of the hydrological cycle based on the period 2003–2012.

An observation-based global land monthly surface air temperature, developed at the Climate Prediction Center, National Centers for Environmental Prediction, is used to calculate the rate of change in the mean annual temperature for the period 1982–2012. The dataset uses station observations collected from the Global Historical Climatology Network version 2 and the Climate Anomaly Monitoring System (GHCN + CAMS)65. The rate of change in temperature is then used to compute the percentage of change in P, ET and Q equivalent to the temperature increase of one degree Kelvin (1°K). For example, the percentage of change in P equivalent to the temperature increase of 1°K can be calculated as \(\frac{{rate\,of\,change\,in\,P}}{{climatology\,P}} \times \frac{{100\% }}{{rate\,of\,change\,in\,temperature}}\).

Global map of drylands and non-drylands based on the first edition of the World Atlas of Desertification66 was used in the analysis of the trends across wet and dry environments. Drylands include ‘Hyper-Arid’, ‘Arid’, ‘Semi-arid’ and ‘Dry Sub-humid’. Non-drylands include ‘Humid’ and ‘Cold’ classes. The classification of the land is based on the seasonal variability of precipitation and an aridity index \({{{\mathrm{AI}}}} = \frac{P}{{PET}}\) where P and PET are the mean annual precipitation and mean annual potential evapotranspiration respectively.

These datasets are listed in Supplementary Table 1.

Precipitation regimes

We cluster the global land into nine precipitation regimes using k-mean unsupervised classification based on the algorithm developed by Lloyd (1957)67. The optimal number of clusters, i.e. nine clusters, is chosen with silhouette analysis68 by computing the average silhouette of predictors over several values of clusters in the range [4–15] and achieving the highest average. The average silhouette approach is commonly used to measure the quality of clustering and for finding the optimal number of clusters69.

We build the k-mean classifier on two predictors. These are mean annual precipitation totals and within-year relative standard deviation (%) both computed and averaged over the period 1980–2018 for MSWEP-P. The derived precipitation regimes map and a scatter plot displaying the pixel values of the two predictors colour-coded by precipitation regime are illustrated in Fig. 1. The global land is clustered into four dry and five wet regimes, and are labelled based on the values of their centroids form the driest to wettest as: very dry with very high variability (V.dry, V.H.variability), dry with high variability (Dry, H.variability), dry with medium variability (Dry, M.variability), mild dry with medium variability (M.dry, M.variability), mild wet with medium variability (M.wet, M.variability), wet with medium variability (Wet, M.variability), wet with low variability (Wet, L.variability), very wet with low variability (V.wet, L.variability), and extremely wet with low variability (Ex.wet, L.variability). The driest and the wettest regimes occupy the largest (32.5%) and the smallest (0.1%) proportion of the land area respectively. We note that using the mean annual apportionment entropy as a predictor instead of the relative within-year standard deviation did not produce noticeable differences in the final precipitation regimes map. In a further step, we use land classification of the World Atlas of Desertification to pull apart ‘Drylands’ and ‘Nondrylands’ regions within the dry precipitation regimes.

Analysis of trends

Annual totals are computed from monthly estimates of the components of the hydrological cycle over 1980–2012. For NDVI, annual averages are computed from monthly estimates over 1982–2012. We use a nonparametric trend test as proposed by Mann-Kendall70,71 and the Sen’s slope method72, and we compute the confidence interval of the slope at 95% significance level. A trend in any climate variable is considered statistically significant if the confidence interval of the slope does not include a mix of positive and negative values. This way, we ensure that statistically significant trends are strictly positive or negative, which is not necessarily guaranteed when considering the p-value. A trend is considered uncertain positive or uncertain negative if its slope is positive or negative respectively, and its confidence interval includes a mix of positive and negative values.

We also identify inconclusive trends using data assimilation, and uncertainties in the individual components of the hydrological cycle:

Data assimilation

We enforce the closure of the hydrological cycle at every grid cell and month during 2003–2012 using a data assimilation technique (DAT)64,73 that modifies each component of the hydrological cycle within its pre-defined uncertainty to achieve the balance. In this analysis, we use the reference ΔTWS is based on GRACE Mascons dataset. Uncertainties in Q, ET and ΔTWS are provided in the data. We use the relative uncertainty of GPCC as indicator of uncertainties in P. We then analyse the trends in each component before and after applying the DAT. If directions of trends in individual components are inconsistent pre- and post-DAT, longer period trends (i.e. 1980–2012) are deemed inconclusive. As a result of this robustness measure, trends in P, ET, and Q have been considered inconclusive at 7%, 5% and 12% of the land respectively.

Uncertainties in the individual components of the hydrological cycle

We incorporate uncertainties in ET and Q to further assess the reliability of the detected trends following these steps for each flux (i.e. ET and Q): (i) We derive 50 different random samples of flux’s monthly time series within the interval flux ± uncertainty; (ii) we derive trends in annual total flux for each sample dataset over 1980–2012; (iii) trends exhibiting consistent direction across less than 90% of the samples are considered inconclusive. Trends in P are compared to trends in the reference precipitation dataset (GPCC), and those showing contradictory change direction across the two precipitation datasets, are considered inconclusive. For NDVI, trends exhibited at low-quality grid cells (i.e. flag = 0) are deemed inconclusive.

Finally, we define ‘robust trends’ positive or negative trends that are statistically significant, and are neither uncertain, nor inconclusive. These measures ensure that our findings are on the basis of a robust trends analysis and that we are not drawing conclusions from spurious trends.