Main

Freshwater systems globally are under unprecedented pressure from human actions. Water extraction and infrastructure, land use and land cover change, and climate change now considerably modify the quantity and timing of atmospheric and terrestrial freshwater flows, with crucial implications for Earth’s climate and ecosystems1. Remarkable signals of water cycle changes include widespread and severe river flow regime alterations2, intensification3 and homogenization4 of the global water cycle, profound changes in terrestrial water storage5, and increases in the severity, frequency and duration of floods and droughts6.

Scientific understanding of the magnitude and nature of anthropogenic water cycle change—and particularly change that is relevant to the hydroclimatic and hydroecological stability of the Earth system—is yet limited. This is due to three main shortcomings. First, water cycle changes are typically analysed over relatively short time periods, often decades7,8,9,10,11, which leads to disregarding natural variability over multidecadal to centennial scales. Second, changes are often quantified against reference conditions that are already affected by anthropogenic influences such as climate change12,13,14,15, hiding the full extent of anthropogenic change. Finally, change is typically assessed in terms of single elements of the water cycle2,8,14,16,17,18,19,20,21 or by aggregating freshwater flows and stocks very broadly (whether in space10,22,23 or in time11,24), which may disregard the interconnectedness of the water cycle and its integral role in the Earth system25,26. Thus, a global assessment that acknowledges the Earth system relevance of freshwater by using long time scales and adequate reference conditions together with an approach incorporating interactions and feedbacks between different water cycle elements is still missing.

In this Article, we quantify deviations of blue (represented by streamflow) and green water (represented by root-zone soil moisture—hereinafter simply soil moisture—that is, the soil water available to plants) from pre-industrial baseline conditions. Jointly assessing change in streamflow and soil moisture conditions, using coherent methodology, allows for qualitatively inferring changes in other water cycle elements and the varying drivers of change. Moreover, our cross-scale assessment method enables analysing change at local, regional and global scales. We establish a baseline state before the onset of major human impacts on the water cycle using the pre-industrial period 1661–1860 as reference conditions and assess change over the following 145 years (1861–2005) against it—illustrating the trajectory of how human-driven change in streamflow and soil moisture has evolved. Widespread deviations from the ‘pristine’ state that are detected by our approach can be considered to pose elevated risks to the Earth system functions of freshwater. For example, terrestrial and freshwater ecosystems have adapted to specific quantities and timing of water flows such that changes in streamflow and soil moisture can severely affect ecosystem status and biodiversity27,28,29. Moreover, the wetness of landscapes regulates climate from micro- to regional and global scales, such that changes in the magnitude and timing of wetness can impact rainfall and consequently streamflow and soil moisture both locally and remotely30.

In addition to a comprehensive standalone assessment of changing streamflow and soil moisture conditions globally, our analysis concludes recent efforts to revise the planetary boundary (PB) for freshwater31,32. PBs set limits to nine processes, including freshwater change, that together regulate the state of the Earth system and thus delimit a safe operating space for retaining a quasi-stable Holocene-like state33,34,35. The recent third assessment of the PB framework uses our approach and globally aggregated results to redefine and quantify a PB for freshwater change33. The new definition replaces the previous estimates of the PB for freshwater use, which have been criticized for their limited capacity to capture the interconnected direct and indirect anthropogenic pressures on the water cycle31. Our approach refines and advances the method first proposed by Wang-Erlandsson et al.32 for quantifying a PB for green water by presenting a coherent, comparable and spatially explicit assessment of blue (streamflow) and green water (soil moisture) change.

The contribution of this paper is, thus, twofold. First, we present a comprehensive, cross-scale assessment of how, where and when streamflow and soil moisture have changed globally from their pre-industrial reference conditions, utilizing an ensemble of consistently forced, state-of-the-art global hydrological models (GHMs; Fig. 1 and Methods). Second, we use our results to unpack the spatial and temporal distribution of streamflow and soil moisture deviations underlying the new freshwater change PB to substantiate and further interpret the globally aggregated results used in the latest PB assessment33.

Fig. 1: Methodological outline.
figure 1

a–f, The analysis steps consist of setting the local baseline range (a), identifying local deviations (b), computing the percentage of land area with local deviations (c), defining pre-industrial variability (d), and comparing the industrial period against the pre-industrial period locally (e) and at global and regional scales (f). The steps a–d are described in detail in Methods.

Full size image

Estimating streamflow and soil moisture deviations

To estimate streamflow and soil moisture change across scales, we used data simulated by an ensemble of state-of-the-art gridded GHMs forced with bias-adjusted CMIP5-generation climate models and dynamic socio-economic conditions (Methods and Supplementary Table 1). The data were given in two time periods: a pre-industrial (1661–1860) and an industrial period (1861–2005).

Based on the 200-year pre-industrial period, we defined the local baseline range of streamflow and soil moisture, separately for each 0.5° grid cell and month, as the range between the 5th (dry bound) and the 95th (wet bound) percentiles of pre-industrial monthly streamflow or soil moisture (Fig. 1a). We then compared pre-industrial streamflow and soil moisture data against the local baseline range to detect local deviations (Fig. 1b). Months with values below the dry bound were marked as dry local deviations and months with values above the wet bound were marked as wet local deviations (Fig. 1b). The areas of grid cells with local streamflow or soil moisture deviations were then summed up globally (or regionally). We divided this sum by the total ice-free land area to yield the percentage of land area with local deviations (Fig. 1c). Finally, to establish global (or regional) reference conditions, we defined pre-industrial variability (Fig. 1d) from the pre-industrial percentage of land area with local deviations. Two key metrics characterize pre-industrial variability: the median and the upper end of pre-industrial variability (Fig. 1d). The median corresponds to the 50th and the upper end to the 95th percentile of the percentage of land area with local streamflow or soil moisture deviations during the pre-industrial period.

After defining the pre-industrial reference conditions, we compared streamflow and soil moisture data during the industrial period against them. We repeated the detection of local deviations (Fig. 1b) and computed the percentage of land area with local deviations (Fig. 1c), which allowed us to compare changes in local deviation frequency (Fig. 1e) and in the percentage of land area with local deviations (Fig. 1f).

Global land area with local deviations

The global land area with local streamflow or soil moisture deviations showed relatively little variation during the pre-industrial period (Fig. 2), as expected because of our modelling setup that uses fixed forcing to prescribe the pre-industrial period as stable (Methods). Therefore, these pre-industrial conditions can be considered a useful reference baseline that existed before the onset of major anthropogenic impacts on the water cycle. The median of pre-industrial variability (that is, the typical land area with local streamflow or soil moisture deviations) was at 9.4% for streamflow and 9.8% for soil moisture (Fig. 2). This is in accordance with the expectation that dry or wet deviations should each occur in 5% of the data points in each grid cell over the pre-industrial period (Methods). Annually, the pre-industrial percentage of land area with local deviations varied mostly within ±1.5 percentage points (p.p.) around the median of pre-industrial variability (Fig. 2).

Fig. 2: Percentage of global ice-free land area with local deviations.
figure 2

a,b, For streamflow (a) and soil moisture (b), the annual percentage is shown, which is computed as an average of monthly percentages. The annotated years mark the 10-year moving (trailing) mean transgressing the upper end of pre-industrial variability (95th percentile of the pre-industrial percentage of land area with local deviations; Fig. 1d). The ensemble median and interquartile range are computed from n = 23 (streamflow) and n = 15 (soil moisture) ensemble members. Values before 1691 are excluded and the ensemble median line for 1861–1890 is shaded and dashed due to traces of model spinups being common during these years (Methods).

Full size image

However, occurrence of streamflow and soil moisture deviations started to steadily increase after the end of the pre-industrial period, with both global indicators surpassing the upper end of pre-industrial variability in the early twentieth century (Fig. 2). Global land area with local streamflow deviations transgressed the upper end of pre-industrial variability (10.2%) already in 1905. The degree of transgression continued to increase, apart from two dips around 1940 and 1970. By the end of our study period (mean of 1996–2005), areas with local streamflow deviations covered 18.2% (~24 million km2) of global ice-free land area (Fig. 2a). This absolute increase of 8.0 p.p. corresponds to a relative increase of 78% compared with the upper end and almost a doubling (8.8 p.p. absolute, 94% relative increase) compared to the median of pre-industrial variability. The evolution of streamflow deviation occurrence is mainly due to dry deviations, which were still increasing in spatial coverage by the end of our study period, while the increase in wet deviations had plateaued decades earlier (Extended Data Fig. 1b,c).

Soil moisture exhibits a similar, though slightly less pronounced, global trajectory and pattern: the upper end of pre-industrial variability (11.1%) was transgressed in 1929, after which the land area with local soil moisture deviations increased consistently, apart from two dips (like streamflow) in the latter half of the twentieth century (Fig. 2b). During the last decade of our analysis period (1996–2005), local soil moisture deviations occurred on 15.8% (~20 million km2) of global ice-free land area, which signifies a 4.7 p.p. absolute (42% relative) increase compared with the upper end of pre-industrial variability and a 6.0 p.p. (61%) increase relative to the median. In contrast to streamflow, by the end of our study period, occurrence of dry soil moisture deviations was declining slowly, yet remaining at considerably higher levels than during the pre-industrial period, while the occurrence of wet deviations continued to increase towards the end of the twentieth century (Extended Data Fig. 1e,f).

When assessing how sensitive our results are to the definition of the local baseline range (Extended Data Fig. 2), we found that transgressions relative to the upper end of global pre-industrial variability become larger when widening the local baseline range. For example, the relative transgression was three- to fourfold larger when using the 1st–99th percentiles as the local baseline range in Extended Data Fig. 2e,f. This indicates that the most extreme conditions increased more than less extreme conditions, but the temporal evolution during the industrial period remains similar regardless of the local baseline definition (Extended Data Fig. 2). Although the model spread is relatively large (Extended Data Figs. 3 and 4), especially towards the present day and in the case of dry streamflow deviations (Extended Data Fig. 1c), the entire ensemble interquartile range of both streamflow and soil moisture surpassed the upper end of pre-industrial variability in the latter half of the twentieth century (Fig. 2). Thus, the global patterns we observed are robust against key methodological sensitivity, and the hydrological models used here mostly agree on the evolution of the percentage of global land area with local deviations. To summarize, there has been—and is still ongoing—a robust and remarkable drift away from the quasi-stable pre-industrial streamflow and soil moisture conditions (characterized by a narrow, ~3-p.p.-wide variability range) to persistently increasing change (currently 4.7–8.0 p.p. beyond pre-industrial variability).

Local and regional evolution of deviations

Mapping changes in the frequency of local deviations reveals a general pattern of more frequent dry deviations of both streamflow and soil moisture in much of the tropics and subtropics, with wet deviations becoming more common in temperate and subpolar regions, and in many highland areas (Fig. 3 and Extended Data Fig. 5). For streamflow, increases in the frequency of dry local deviations are more common, while for soil moisture, neither dry nor wet deviation frequency increases dominate. Comparing local deviation frequency in 1976–2005 against 1691–1860, 39.2% of land area shows a statistically significant increase in dry streamflow deviation frequency, while 9.2% of land area exhibits a statistically significant increase in wet streamflow deviation frequency. For soil moisture, 26.9% of land area shows an increase in dry and 19.2% an increase in wet deviation frequency. Areas in which the frequency of both dry and wet local deviations increased significantly cover 2.0% (streamflow) and 1.0% (soil moisture) of land area. Therefore, 46.4% of land area has experienced a statistically significant increase in streamflow deviation frequency and 45.1% in soil moisture deviation frequency. Considering that 5% of all months were marked as local deviations during the pre-industrial period (Methods), increases in local deviation frequency were major (more than a 5 p.p. increase, that is, more than a doubling) in approximately half of the area with increased frequency of dry streamflow deviations, a third of the area with increased frequency of soil moisture deviations (both dry and wet), and a quarter of the area with increased frequency of wet streamflow deviations (Fig. 3). Increases in the frequency of wet local deviations are concentrated in the boreal winter months (DJF) while the boreal summer months (JJA) are characterized by increased dry local deviation frequency, especially around the Mediterranean, the Sahel region and Central and North America (Extended Data Fig. 6).

Fig. 3: Statistically significant increases in dry and wet local deviation frequency.
figure 3

a,b, For streamflow (a) and soil moisture (b), changes in the frequency of local deviations are computed by comparing the ensemble median frequency of local deviations during 1976–2005 against 1691–1860, and significance of change is tested at a 95% confidence level (P = 0.05) with R package stats function prop.test96 (one-sided test). Colours denoted with ‘+‘ indicate statistically significant increases with magnitude ≤5 p.p. (minor), whereas colours denoted with ‘++’ indicate statistically significant increases with magnitude >5 p.p. (major). Colours denoted with ‘+(+)’ pool together any statistically significant increase (minor or major).

Full size image

Persistent streamflow and soil moisture changes and their timing are spread out unevenly (see Fig. 4, in which local streamflow and soil moisture deviations are aggregated at river basin scales instead of the global scale). In many regions—especially around mid-latitudes—land area with local deviations persistently transgressed the region-specific upper end of pre-industrial variability already before 1940, while in other areas—particularly in humid regions—a transgression has not occurred yet. The Mississippi, Indus and Nile basins, for instance, were among the first regions where persistent transgression happened in the case of streamflow (in the Nile, also for soil moisture) (Fig. 4a,b). For soil moisture, a persistent transgression has occurred in fewer regions and often later than in the case of streamflow (Fig. 4c,d). This is visible particularly in Siberia, South/Southeast Asia and the Congo Basin, where the persistent transgression year for soil moisture followed 10–20 years after that of streamflow—or the transgression has not occurred yet (Fig. 4).

Fig. 4: Timing of persistent transgression of the region-specific upper end of pre-industrial variability.
figure 4

a–d, For streamflow (at HydroBASINS level 2, a, and at level 3, b) and soil moisture (at HydroBASINS level 2, c, and at level 3, d), local deviations are first aggregated regionally to get pre-industrial percentage of land area with local deviations (Fig. 1c), then the upper end of pre-industrial variability is set to the 95th percentile of the percentage of land area with local deviations (Fig. 1d), and finally the percentage of land area with local deviations is calculated for the industrial period and compared against pre-industrial variability (Fig. 1f). The persistent transgression year is defined as the first year, during which the 10-year moving (trailing) mean of percentage of land area with local deviations has exceeded the region-specific upper end of pre-industrial variability for ten consecutive years, without returning below this limit after this year. The regions shown here depict basins delineated by the HydroBASINS data set97 level 2 (a and c; n = 60, mean area 2,247,000 km2, median area 2,045,000 km2) and level 3 (b and d; n = 265, mean area 509,000 km2, median area 312,000 km2).

Full size image

Interpreting local deviation patterns

Assessing the changes in local streamflow and soil moisture deviations in unison and qualitatively comparing them with literature provides more detailed insights into their potential causes and impacts. These qualitative assessments should be considered preliminary, and analytical driver attribution with improved data is warranted in future studies (Methods). The general pattern of increasing wet deviation frequency of both streamflow and soil moisture in the Northern Hemisphere probably arises from the precipitation change following ~1 °C of mean global warming6, which affects both streamflow and soil moisture similarly if major human impacts on land surface, such as land cover change and soil degradation, do not change precipitation partitioning. Another distinctive example of climate-induced deviations is found in the Sahel region (stretching across the African continent south of the Sahara; partially covered by the Niger River Basin in Fig. 5). Drier conditions (relative to pre-industrial climate) dominated in the region already in the first half of the twentieth century36—which agrees with our analysis showing increased dry deviation frequency already when comparing the years 1931–1960 against 1691–1860 (Extended Data Fig. 7)—and intensified from the 1970s onwards. Recent decades in the Sahel have additionally witnessed widespread wet local deviations (and different combinations of dry and wet), which could be explained by changes in temporal dynamics, such as intensity of precipitation, number of wet days, occurrence of dry spells and timing and length of the rainy season36,37,38.

Fig. 5: Combined increases in local deviation frequency, and percentages of land area with local deviations in selected regions.
figure 5

The classification is based on local streamflow and soil moisture deviation frequency increases shown in Fig. 3, and it pools together both minor and major increases (represented by ‘+(+)’ in the legend). Computing the percentage of land area with local deviations (Fig. 1c) was performed within each region, and the median and upper end of pre-industrial variability (Fig. 1d) were also computed regionally. Selected region outlines are derived from the HydroBASINS data set97.

Full size image

As an example of direct human drivers, increased frequency of dry local streamflow deviations coinciding with increased wet local soil moisture deviations in a given area probably indicates an effect of irrigation expansion. Indeed, by inspecting local deviation frequency on areas equipped for irrigation, we found that irrigation intensity is associated with increasing dry streamflow and wet soil moisture deviation frequency (Supplementary Fig. 1). These cases are observed in many heavily irrigated regions39, such as South Asia, eastern China, Western United States and the Nile Delta (Fig. 5, pink colour). These regions are also the likely early drivers of globally aggregated deviations, as their irrigation extent was relatively large already in the early twentieth century39, and local deviations were prevalent already then (Extended Data Fig. 7). Increases in the frequency of both wet and dry streamflow deviations, as well as wet soil moisture deviations, (Fig. 5, dark blue) could indicate yet more signals of river flow regulation (decreased flood peaks and increased dry season streamflow) in combination with extensive irrigation (decreased streamflow and increased soil moisture)23. These cases can be found in many heavily modified river basins40, such as the Nile, the Aral Sea and parts of India and Thailand (Fig. 5).

Our findings of increased streamflow and soil moisture deviation frequency are consistent with many regional to global freshwater-mediated impacts that have been reported independently. One of the most dramatic examples of blue water change is the Aral Sea, where the overuse of water for irrigation resulted in considerable lake depletion, consequent ecological degradation and regional climate change41. In our results, the region’s two large river basins, Amu Darya and Syr Darya, exhibit a steep increase in dry streamflow and a moderate increase in wet soil moisture deviations, starting in the late 1960s (Fig. 5)—in accordance with the substantial irrigation expansion that started in the region at that time41.

Soil moisture changes have been associated with productivity loss, as exemplified by drying-induced forest dieback42 in many of the regions where we find a major frequency increase in dry soil moisture deviations, such as the Mediterranean basin, central North America and West Africa (for example, Niger River Basin in Fig. 5). Productivity shocks due to dry and wet soil moisture deviations (droughts and floods) have also been reported on many cultivated lands, particularly in South and East Asia, Australia and North Africa43, where we find both drying and wetting (Fig. 3b, India in Fig. 5). Other examples of ecological and climatic impacts of water surpluses are habitat loss in the Central Amazon floodplains due to anthropogenic flood pulse disturbance by dams44, and increased greenhouse gas emissions from reservoirs45, though these freshwater alterations are relatively small scale and thus not clearly visible in our maps.

The new planetary boundary for freshwater change

Our approach has been adopted as the new standard for defining and assessing the PB for freshwater change33, replacing the relatively narrow outlook on global freshwater change achieved by previous freshwater use PB estimates. The previous approaches—which assessed water use and withdrawals—fell short in capturing freshwater change due to, for example, land use and land cover change and anthropogenic greenhouse gas and aerosol emissions, which alter evaporation, soil moisture, precipitation and runoff patterns31,46,47,48,49,50. They have also been criticized for aggregating freshwater fluxes too broadly, for example, by operating on an annual scale only or simply summing global water use and availability despite diverse local impacts51,52. Most recently, safe and just Earth System Boundaries were proposed for several domains, including freshwater15—partly building on earlier iterations of the PB framework. However, accommodating for the additional justice aspect in the freshwater Earth System Boundaries resulted in conceptual and methodological trade-offs similar to those of the earlier freshwater use PB approaches (Supplementary Text). The approach presented in this paper captures the breadth of global water cycle change better than the earlier approaches by analysing two key components of the freshwater cycle, by using metrics that relate to freshwater’s Earth system impacts, and by adopting meaningful baselines and long time scales. Although streamflow and soil moisture do not explicitly represent all aspects of the freshwater cycle, such as groundwater, they are connected, either directly or indirectly, to major anthropogenic modifications of the freshwater cycle and associated Earth system functions (Methods and Supplementary Text).

The upper end (95th percentile) of pre-industrial variability in global land area with local streamflow (blue water) and soil moisture (green water) deviations are now the two components of the PB for freshwater change33. This definition represents global conditions occurring expectedly once in 20 years in the pre-industrial reference state and corresponds to 10.2% of the global land area for blue water (Fig. 2a) and 11.1% for green water (Fig. 2b). While these limits should not be considered global tipping points33, their persistent exceedance marks a departure from known, safe pre-industrial conditions (representing longer-term Holocene-like conditions; Supplementary Text), which we consider to pose elevated risks to freshwater’s Earth system functions. According to our results, the freshwater change PB has been transgressed for a long time—the current status is 18.2% of global land area for the blue and 15.8% for the green water component (Fig. 2). This is in stark contrast with the first and second PB framework assessments, which considered the freshwater use PB’s status safe34,35. Our current assessment, however, agrees with other, more nuanced studies (that better consider environmental flows and green water, for example) on the freshwater use PB23,53,54, which considered the boundary to be transgressed.

It should be noted that major uncertainties are associated with setting the freshwater change PB (Methods and Supplementary Text). While data-related uncertainties have been mitigated by using a large ensemble of state-of-the-art GHMs that have been validated against observations55,56,57,58,59, perhaps the most important uncertainty stems from the lack of robust quantifications of Earth system-wide responses to water cycle modifications that surpass pre-industrial variability. These responses may be too complex to be quantified with our current knowledge, and it is also possible that the level of freshwater modifications that drives state changes and increases risks in the Earth system can only be determined in retrospect.

Nevertheless, because the freshwater cycle upkeeps many life-supporting Earth system functions, risks of environmental degradation can be assumed to elevate along with change away from a stable, Holocene-like state1—a change that has occurred in streamflow and soil moisture conditions within the last century (Fig. 2). While it is possible that buffers or stabilizing feedbacks in the Earth system (such as CO2 fertilization) allow for increased variability in the water cycle without increasing risks, there is also evidence that self-reinforcing feedbacks (such as moisture recycling60 and methane61 feedbacks) may amplify the impacts of freshwater change on Earth system functioning. Thus, considering our limited capabilities to model complex, interacting Earth system processes (and even the different aspects of the water cycle), and freshwater’s central role in mediating Earth system interactions25,26, the precautionary principle of the PB framework33,34,35 motivates a conservative placement of the new PB for freshwater change at the lower level of scientific uncertainty. Here and in the recent third PB assessment33, this is considered to be the upper end of pre-industrial variability. The major transgression and remaining uncertainties of this boundary placement, however, warrant continued research on the role of freshwater in the Earth system.

Concluding remarks

We have presented here the trajectory of anthropogenically driven change in streamflow and soil moisture since the pre-industrial period and shown that this change has been pervasive across spatial and temporal scales. Globally, land area in which streamflow and soil moisture deviates from local pre-industrial reference conditions has increased by 78–94% and 42–61%, respectively, within the 145-year industrial period. The patterns we observe in these two key freshwater cycle elements are in accordance with the few other studies reporting long-term (>100 years) changes in other water cycle elements, such as groundwater depletion contributing to 4–9% of total sea-level rise18,21,62, and at least a 21% loss of wetlands since the year 1700 (refs. 17,19). Our results align well also with studies reporting decadal changes and trends in increasing hydrological extremes6,8. Evidence of freshwater-triggered ecological and climatic shifts has been mounting while key components of the freshwater cycle have moved further away from their pre-industrial variability at an alarming rate63,64,65. Our findings indicate a transgression of the PB for freshwater change already around the mid-twentieth century, while climate change66, deforestation67 and many other human pressures on the water cycle continue to pose a major risk of further change. Decreasing these pressures by, for example, committing to ambitious climate action, halting deforestation and respecting environmental flows in water use and management is thus imperative to safeguard the life-supporting functions of freshwater.

Methods

Data selection

Our main data source was the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) data repository (available at https://data.isimip.org, last accessed 30 March 2023), from which we used data of the ISIMIP 2b simulation round experiments68 (Supplementary Table 1). We used root-zone soil moisture (hereinafter soil moisture; ISIMIP output variable ‘rootmoist’) to represent green water, and river discharge (hereinafter streamflow; ISIMIP output variable ‘dis’) to represent blue water as two key components of the freshwater cycle (for further elaboration on variable selection, see Supplementary Text). The models in ISIMIP 2b have been validated against observations, showing adequate performance especially when estimates from individual models are combined in an ensemble modelling approach55,56,57,58,59. By constructing ensembles as large as possible, given the availability of different simulation scenarios, we can tackle the uncertainty to which global modelling is always subject to, to a certain degree69.

Two simulation scenarios were required to first determine the pre-industrial reference state (Fig. 1a–d) and then to compare the industrial streamflow and soil moisture conditions against it (Fig. 1e,f). For the pre-industrial period (1661–1860), we used model outputs forced with ‘picontrol’ climate and ‘1860soc’ land use and socio-economic conditions, which described our pre-industrial reference state with a constant 286 ppm CO2 concentration and fixed pre-industrial land use and socio-economic conditions68. The pre-industrial simulation years are, however, nominal, as the CMIP5 data used to force ISIMIP 2b models represent 200 control years with climate variability at pre-industrial levels (with no correspondence to the actual weather of individual years)68. Therefore, while the pre-industrial simulations can be considered to describe an artificially stable state, they adequately correspond to long-term pre-industrial conditions before major anthropogenic modifications affecting the distribution of freshwater. At the transition from pre-industrial to industrial time (as defined by the ISIMIP protocol) in 1860, we switched to ‘historical’ climate combined with ‘histsoc’ land use and socio-economic conditions, in which carbon-induced climate change and human influences (for example, land use change, water use and dam operation) were dynamically represented in the data68. The ISIMIP 2b industrial period covers years from 1861 to 2005, the end year owing to estimates of, for example, irrigation extent not being available39.

In ISIMIP 2b, streamflow and soil moisture were simulated by GHMs, which were forced with bias-adjusted output of modelled climate from general circulation models (GCMs). The ISIMIP 2b modelling groups have bias-adjusted the GCM outputs from CMIP5 models70 with EWEMBI, a data product combining climate re-analysis and observations, covering the years 1979–2013 (ref. 68). We selected all GHM–GCM combinations, for which output data were available for both picontrol–1860soc and historical–histsoc scenarios (Supplementary Table 1). For soil moisture, outputs from four GHMs (CLM50 (ref. 71), LPJmL72, MPI-HM73 and PCR-GLOBWB74) were available, and for streamflow, outputs from six GHMs (H08 (ref. 75), LPJmL, MATSIRO76, MPI-HM, PCR-GLOBWB and WaterGAP2 (ref. 77)) were available. For all GHMs, outputs of four GCMs were available: GFDL-ESM2M78, HadGEM2-ES79, IPSL-CM5A-LR80 and MIROC5 (ref. 81), except for MPI-HM, which lacked data forced with HadGEM2-ES for both streamflow and soil moisture. Hence, our soil moisture data ensemble consisted of 15 members (GHM–GCM combinations) and our streamflow data ensemble consisted of 23 members. For streamflow, we aggregated daily values to monthly values by taking the mean of daily streamflow, whereas soil moisture data were readily delivered in monthly time resolution, which is the highest temporal resolution available in ISIMIP 2b for soil moisture. The spatial resolution of all data was 0.5 degrees (ca. 50 × 50 km at the Equator). These temporal and spatial resolutions are standard practice in global hydrological modelling, partially stemming from model input data resolutions68,82, while the development of hyper-resolution global hydrological modelling is in its early stages (see, for example, ref. 83).

Data preparation and intercalibration between periods

Because the ISIMIP scenario simulations can be run independently of each other, the end state of the pre-industrial simulation does not necessarily equal the initial state of the industrial simulation at the grid cell scale. Two kinds of discontinuities may thus result: (1) traces of model spinup, that is, GHMs reach a water storage equilibrium with a delay after the beginning of the simulation, and (2) abrupt shifts in value distribution around the transition point at the end of year 1860. Therefore, we performed an intercalibration between the pre-industrial and industrial periods with a two-step approach applied separately for each grid cell, month and ensemble member. Supplementary Fig. 2 shows an example of traces of model spinup, whereas Supplementary Fig. 3 shows our intercalibration approach along with an example of a cell in which the value distribution shifts. This intercalibration is additional to and separate from the GCM bias adjustment conducted by the ISIMIP modelling groups68.

First, we checked for traces of spinups by detecting the most likely changepoint in the mean and variance with function cpt.meanvar in R package changepoint84,85 separately for the pre-industrial and industrial periods. As the trace of spinup relates to reaching a hydrological equilibrium in a cell, we assumed that a true trace of spinup is located between 10 and 30 years after the initiation of each simulation period. Outside this period, we considered that detected changepoints were true (natural) changes in mean and variance and did not imply a trace of spinup. If a changepoint was detected between 10 and 30 years after initiation (as in Supplementary Fig. 2, for example), all values before the changepoint were excluded from the intercalibration; otherwise (as in Supplementary Fig. 3, for example), values from the ten first years of the simulation periods were discarded. Depending on the month and ensemble member, changepoints indicating traces of spinup were detected in 1.4–8.7% (pre-industrial streamflow), 1.6–7.2% (industrial streamflow), 0.4–4.6% (pre-industrial soil moisture) and 0.6–5.7% (industrial soil moisture) of all cells.

Second, we checked and corrected shifts in value distribution around the simulations’ transition point in 1860, using an iterative technique combining linear extrapolation and quantile mapping. Because the simulation periods have no temporal overlap, we extrapolated pre-industrial data onto the industrial simulation period (that is, extended the pre-industrial time series to nominally cover years past 1860). This was done by fitting a linear regression model86 to pre-industrial values excluding spinup and extrapolating the linear trend (blue lines in Supplementary Fig. 3a). Further, we computed the standard deviation of the pre-industrial values to which the linear model was fit (defined here as σpreind) and added normally distributed random noise (μ = 0, σ = σpreind) around the linear extrapolation line to create extrapolated data points (blue circles in Supplementary Fig. 3a).

After extrapolation, we fitted non-parametric quantile mapping using robust empirical quantiles implemented in R package qmap87,88 between industrial simulation values and extrapolated pre-industrial values. Four quantiles were used in fitting the quantile mapping. We treated the industrial simulation values as ‘observed’ data and the extrapolated pre-industrial values as ‘modelled’ data. We chose to fit quantile mapping between 30 years of industrial and extrapolated pre-industrial values, starting from the first year not excluded after spinup detection (pink shading in Supplementary Fig. 3a). Hence, quantile mapping was mostly fitted using years 1871–1900 as the ‘observed’ values but with flexibility of extending the fitting period to 1891–1920 at maximum. For example, the fitting period for the cell shown in Supplementary Fig. 2 was 1874–1903, whereas for the cell shown in Supplementary Fig. 3, the fitting period was 1871–1900. Although 30 years is a relatively short period, and some spinups might end only after 1891, we chose to limit the quantile mapping with these years to prevent the industrial values being substantially affected by anthropogenic impacts, which increases in likelihood during the twentieth century.

Finally, pre-industrial values excluding spinup were corrected with the fitted quantile mapping function (Supplementary Fig. 3b). However, because distinctive individual values or undetected traces of spinup may have distorted the quantile mapping function, we checked whether applying quantile mapping succeeded in improving the fit between the simulation periods. To do this, we ran our analysis up until detecting local deviations (Fig. 1a,b)—first without correction and then with correction—and checked whether the number of local deviations (Fig. 1b) decreased owing to the correction. Detection of local deviations was performed for a period of 50 years, beginning from the first year not excluded after spinup detection (hatched fill in Supplementary Fig. 3). This was to include data outside the quantile mapping fitting period. If the number of local deviations decreased, we considered that our intercalibration had improved the fit between the simulation periods and the corrected cell values were accepted; otherwise we retained the uncorrected data. Globally, depending on month and ensemble member, the correction was accepted for 6–26% of streamflow cells and 10–29% of soil moisture cells. It should be noted that no constraints were set to preserve the lateral distribution of water in the intercalibration, but as soil moisture or streamflow changes are expressed in our analysis only relative to values in the same grid cell (that is, we do not compare values across cells, nor are we routing water), this has a negligible impact on our key results. For globally aggregated results, the intercalibration has the most impact in the beginning of the industrial period, by decreasing the global land area with local deviations, while towards the end of our time series, the results with and without intercalibration are nearly in agreement (Supplementary Figs. 4 and 5).

Setting the local baseline range and identifying local deviations

In setting the local baseline range (Fig. 1a), identifying local deviations (Fig. 1b) and computing the percentage of land area with local deviations (Fig. 1c), we followed the general approach of Wang-Erlandsson et al.32 with some modifications. Like them, we set the 5th and 95th percentiles of pre-industrial streamflow and soil moisture to bound the local baseline range (separately for each grid cell, month and ensemble member), but we drew these bounds from the considered grid cell’s values only, whereas Wang-Erlandsson et al.32 drew their local baseline range from the values in the considered grid cell and its neighbourhood. We chose a stricter definition of the local baseline range because including neighbourhood values in determining the local baseline range can potentially set it unrealistically wide—especially in the case of streamflow, in which neighbourhood values strongly depend on flow directions.

Spinups were excluded from setting the local baseline range, which means that for most cells, dry and wet bounds (Fig. 1a) were drawn from 190 values covering years 1671–1860, though allowing flexibility up until 1691 (for example, for the grid cell shown in Supplementary Fig. 2, the bounds were drawn from values in 1690–1860). We also evaluated the sensitivity of our results to the definition of the local baseline range by repeating our main global analysis (Fig. 2) using the 2.5th and 97.5th and the 1st and 99th percentiles as the dry and wet bounds (Extended Data Fig. 2).

Then, again for each grid cell, ensemble member and pre-industrial month, we checked whether streamflow or soil moisture fell within the local baseline range, to identify local deviations (Fig. 1b). In case a monthly value was lower than the dry bound (5th percentile of pre-industrial values), this month was marked as a dry local deviation. Contrarily, in cases of the monthly value surpassing the wet bound (95th percentile of pre-industrial values), the month was marked as a wet local deviation. Local deviations were evaluated for all months including spinups, and all deviations were determined in a binary fashion, as ‘deviated’ or ‘not deviated’, with no regard to deviation magnitude.

Computing the percentage of land area with local deviations and defining pre-industrial variability

After identifying local deviations (Fig. 1b), we summed up the land areas of grid cells with dry and wet local deviations, separately for each month and ensemble member. We excluded Antarctica and other permanent land ice areas using HYDE 3.2.1 anthromes89 and divided the land area with local deviations by total ice-free land area to yield the percentage of land area with local deviations (Fig. 1c). We transformed monthly percentages into annual percentages by taking annual means. This differs from Wang-Erlandsson et al.32 who marked an individual grid cell as deviating during a given year if a local deviation occurred in any month of that year, and computed the global percentage of land area with local deviations only annually. We chose to do this at a monthly timestep to enable seasonal analysis (Extended Data Fig. 6) and to distinguish between cases in which local deviations occur in only one month or many months within a year. Thus, local deviations that spread out over multiple months of the year have a higher impact on the annual percentage of land area with local deviations.

As approximately 5% of all pre-industrial values in each grid cell were marked as dry and wet local deviations by the definition of the local baseline range (5th–95th percentiles), the expected global percentage of land area with local deviations was also approximately 5%, for dry and wet deviations separately (under conditions with little interannual variance, that is, not during spinups). Hence, summing the percentages of land area with dry and wet local deviations together, the expected value for the global percentage of land area with local deviations was approximately 10%. As spinups were prevalent especially in some GHMs (Extended Data Figs. 3 and 4), global percentages of land area with local deviations before 1691 were excluded from further analysis. Therefore, we ended up with a time series of n = 170 for the annual pre-industrial percentage of land area with local deviations.

The pre-industrial percentage of land area with local deviations (Fig. 1c) was further used to define pre-industrial variability (Fig. 1d). We took the ensemble median (n = 15 or n = 23) of the percentage of land area with local deviations to define two main metrics: the median of pre-industrial variability and the upper end of pre-industrial variability (Fig. 1d). The median of pre-industrial variability was defined as the 50th percentile of the percentage of land area with local deviations, while the upper end of pre-industrial variability was defined as the 95th percentile of the percentage of land area with local deviations. Hence, the median describes the typical percentage of land area with local deviations in the pre-industrial (presumably Holocene-like; see Supplementary Text) reference state, whereas the upper end describes conditions that occurred expectedly once every 20 years.

For the industrial period (1861–2005), we identified local deviations (Fig. 1b), computed the percentage of land area with local deviations (Fig. 1c) and compared that with pre-industrial variability (Fig. 1d). We included all years in identifying local deviations but chose to shade years 1861–1890 in the presented results, owing to spinups being prevalent in some GHMs (Extended Data Figs. 3 and 4). Finally, we took the ensemble median of the global percentage of land area with local deviations and applied a moving (trailing) 10-year mean over this ensemble median time series, with the latest 10-year mean (for year 2005, computed as a mean of values from 1996–2005) being defined as the current status.

Limitations and uncertainties

GHMs are known to be sensitive to their implementation of anthropogenic drivers and impacts57,69, of which especially water use, including irrigation, directly affects streamflow and soil moisture but with considerable uncertainty90. ISIMIP 2b simulations attempt to minimize uncertainties stemming from this by using state-of-the-art input data and consistent scenario definitions68. Nevertheless, the different process descriptions of GHMs lead to a considerable spread in our results (Extended Data Figs. 3 and 4), which means that our global—and especially local—results are subject to noteworthy uncertainty. However, we chose not to perform an explicit validation for our hydrological data because ISIMIP 2b data have shown adequate55,56,57,58,59 (though variable91) performance against observations, observed data to validate pre-industrial values are unavailable, and the ISIMIP 2b hydrological data used here are forced by GCM outputs instead of observed climate.

Due to multiple compounding uncertainties, we use a maximally large data ensemble to capture the range of possible outcomes. We report the ensemble median as our best estimate as more sophisticated ensemble construction methods are unavailable in the absence of observations and the ensemble median is often an adequate choice55,58. Moreover, our estimate of the ‘current’ status ends in 2005, due to the ISIMIP 2b simulation protocol68. Since the end of our study period, global trends in many of the key drivers of freshwater change, such as irrigated area92, water use93, dam construction94 and forest loss95, have increased, and therefore, the results presented here are probably conservative.

As our scenario setup of using dynamic histsoc socio-economic conditions against static 1860soc (Supplementary Table 1) already implies a change in anthropogenic drivers of water cycle change, it is not an unexpected result that aggregate streamflow and soil moisture changes were manifested in the early twentieth century (Figs. 2 and 4, and Extended Data Fig. 7). Repeating the analysis using data that are absent of anthropogenic forcing (ISIMIP 'nosoc' scenarios68) would potentially aid in estimating how large proportions of the change shown here are due to direct (for example, water use) or indirect (for example, climate change) anthropogenic factors. However, model outputs of both picontrol–1860soc (or picontrol–nosoc) and historical–nosoc scenarios, using the same GHMs and GCM forcing, were not available in ISIMIP 2b. Therefore, analytical attribution of drivers of the detected change is left for future studies. Future modelling efforts, including the ISIMIP3a simulation round82 (which is forced with more advanced CMIP6 climate data), should provide an excellent starting point for these attribution analyses.

Despite the related uncertainties, our approach is as coherent as current knowledge and modelling capacities allow for a global study. Moreover, the magnitude and rate of change in streamflow and soil moisture conditions in the industrial period and the presumably conservative estimation of the freshwater change PB’s current status suggest that our conclusion of streamflow and soil moisture substantially transgressing their pre-industrial variability ranges is valid.