Most River Basins will Follow their Budyko Curves under Global Warming

,


Introduction
Assessing future shifts in water resources and secure these resources through adaptation and mitigation requires an understanding of hydroclimatic change (Baldassarre et al., 2015;Brown et al., 2019;Nissan et al., 2019;Sivapalan and Blöschl, 2015).For decades, the Budyko framework (Budyko, 1974(Budyko, , 1948) ) has been used to understand hydroclimatic change by studying the relationship between water and energy available on the land surface and considering evaporation's water and energy limits.The framework provides curvilinear relationships for a given hydrological basin, known as Budyko-type curves, between the evaporative ratio (i.e., actual evaporation over precipitation) and the aridity index (i.e., potential evaporation over precipitation) at mean annual or longer scales.Every basin on Earth has a set of combinations of evaporative ratio and aridity index related to its vegetation, soils, topography, climate seasonality, and snow to rain ratio, conforming a Budyko curve.
The change in the location in Budyko space of a particular region or basin can be referred to as "movement in Budyko space" and corresponds to the joint change in the aridity index and evaporative ratio between two time periods (Destouni et al., 2013;Jaramillo et al., 2018;Jaramillo and Destouni, 2014;van der Velde et al., 2013).The movement in the Budyko space can help estimate the effects of particular drivers of change in the evaporative ratio.For example, under stable landcover and water storage conditions, a hydrological basin is expected to move along its Budyko curve if mean-annual changes in the aridity index are the main driver of changes in the evaporative ratio.However, in a more common scenario, other drivers than aridity, such as land cover changes, water use, changes in water storage or snow to rain ratios, or even in the seasonality of both precipitation and potential evaporation may change the evaporative ratio (Jaramillo and Destouni, 2014).For instance, Donohue et al. (2007) state that under stable conditions of the aridity index, land conversion from forest to grassland decreases evaporation and root zone capacity, eventually increasing runoff and decreasing the evaporative ratio (Nijzink et al., 2016;Sterling et al., 2013).In other words, a basin experiencing such widescale land conversion would move vertically downwards in the Budyko space, in the absence of any other driver.On the contrary, land conversion from grassland to forest cover would increase evaporation, would typically move a basin vertically upwards in the Budyko space.Other examples of upward movements in Budyko space include the expansions of irrigation (Wang and Hejazi, 2011) or the impounding effects of reservoirs on rivers (Levi et al., 2015).Furthermore, the non-stationarity conditions of energy and water availability reflected by changes in the snow to rain ratios, shifts in the precipitation regime, or the seasonality of energy availability in a given basin may also drive a movement in Budyko space.For example, a change from snow towards rain decreases runoff (Berghuijs et al., 2014), also translating into an upward movement in the Budyko space.
Hence, the combined effect of changes in the aridity index and any of these different drivers may move basins in Budyko space along trajectories that deviate from the Budyko curves.For instance, on a global scale, Jaramillo and Destouni (2014) studied almost 900 river basins worldwide during 1901-2018 to find that 74% deviated from their potential Budyko curves, that is, they moved beyond the range of 45 to 90 degrees that charactirezies the slope of any Budyko curve (Fig. 1).This result pointed to the additional effects of other drivers such as land and water use, changes in water storage or precipitation seasonality in the river basins on the evaporative ratio.
The question to what extent these others drivers of change (besides the aridity index) will be still significant with ongoing and accelerating global warming.For example, can we expect future movements of hydrological basins under global warming to occur along their Budyko curves, and how will these compare to historical movements?The question becomes increaslingly relevant with the growing effects of greenhouse gas emissions on Earth's water cycle (Gudmundsson et al., 2021(Gudmundsson et al., , 2017;;Huntington, 2006) and its intensification (Huntington, 2006;Koutsoyiannis, 2020).
Here we test the hypothesis that most river basins in the world will move along their Budyko curves due to a dominant effect from increasing aridity index as potential evaporation increases due to global warming.We first calculate movements in Budyko space under historical and future hydroclimatic change from 1900 to 2100 for 353 large river basins worldwide based on hydroclimatic outputs from nine models from the Coupled Model Intercomparison Project -Phase 5 (CMIP5).We then compare these simulated Budyko-space movements with those expected from long-term theoretical and analytically-derived changes in the aridity index.Finally, we discuss the causes and implications of our findings.

Hydroclimatic data
We selected nine Earth System Models that provide monthly data at a high spatial resolution on net radiative forcing (Rn; surface downwelling shortwave radiation), precipitation (P), and actual evaporation (E) (Table 1).These data are required to estimate potential evaporation (E0), the aridity index (E0/P) and the evaporative ratio (E/P) for any location on the Earth's land surface.This selection of models has been previously used to assess hydroclimatic change over the African continent (Piemontese et al., 2019) and excludes models with resolutions coarser than 2.5 degrees to reduce the use of models with poor performance, especially in coastal hydrological basins.Furthermore, all CMIP5 models use the same present and future land cover and land-use scenarios based on Hurtt et al. (2011).
The CMIP5 data was downloaded from the Earth System Grid (ESG, http://www.earthsystemgrid.org),using the realization r1i1p1 (r: realization, i: initialization, p: perturbation) as it provides the largest number of simulations (Taylor et al., 2012).In addition, E0 was estimated using the original energy-only method as it is the best suited for climate model outputs (Milly and Dunne, 2016) (Eq.1).The E0 estimate is expressed as where Rn is the water equivalent of net radiation (net radiation divided by the latent heat of vaporization) in mm/day and G the heat flux into the subsurface (also in mmd/day).We used the surface downwelling shortwave radiation output of the models (rsds) as the net radiation and assumed G to be zero since mean annual potential evaporation is mainly insensitive to seasonal variations (Cook et al., 2014).The constant 0.8 reflects the fraction of available energy (~80%) going into latent heat flux (Koster and Mahanama, 2012).The estimates of P and E estimates are the outputs of the CMIP5 models, pr and evpbs, respectively, and coastal grid cells were eliminated to remove the effect of ocean evaporation on evpbs. of the movement from period 1 to period 2 is the change in the aridity index (E0/P), Δ(E0/P), while the vertical component is the change in the evaporative ratio (E/P), Δ(E/P).The hypotenuse of both components yields the intensity of movement (I), and the angle () represents the direction of movement.The direction (b) and magnitude (Ib) of movement along the Budyko curve represents long-term changes in the aridity index (Eq.4; Yang et al., 2008).Dashed lines represent the energy and water limits constraining hydroclimatic conditions of energy and water availability.The 353 hydrological basins selected for the study are the largest available in the Global Runoff Database Centre GRDC (grdc@bafg.de).All climatic variables were calculated based on the average of all the pixels within each river basin.

Movement in Budyko space from CMIP5 projections
We quantified terrestrial hydroclimatic change by defining four 30-year periods: 1910-1939, 1961-1990, 2010-2039 and 2070-2099.We defined hydroclimate changes as the magnitude and direction of movement in Budyko space between any of these periods, know referred to "change periods", caused by changes in E0/P and E/P.We determined the aridity index and evaporative ratio based on the 30-year averages of E0, E and P, calculated or obtained directly from the projections of the historical experiment and two of CMIP5's Representative Greenhouse Gas Concentration Pathways (RCP4.5 and RCP8.5).The historical experiment is based on simulations forced by observations of atmospheric conditions and accounts for land cover manuscript submitted to Water Resources Research changes during the twentieth century (Moss et al., 2010).RCP4.5 corresponds to the midrange mitigation emission scenario adopted by the Paris Agreement, while RCP8.5 reflects the highest emission rate scenario without carbon emission mitigation strategies.We focused on RCP8.5 since it results in the largest imbalance in the Earth's radiative budget by 2100 (i.e., 8.5Wm −2 at the top of the atmosphere) and therefore represents the largest impacts of carbon emissions and global warming on the water cycle.
We determined the direction and magnitude of the vectors representing movement in Budyko space for the change periods 1910-1939 to 1961-1990, 1961-1990 to 2010-2039, and 2010-2039 to 2070-2099.The movement vector ( ⃗) represents the hydroclimatic change experienced by any hydrological basin over time, with direction (θ), magnitude (I), and horizontal and vertical components of change in E0/P and E/P, Δ(E0/P) and Δ(E/P), respectively.
Some applications of the Budyko framework plot the runoff coefficient (R/P) instead of E/P on the y-axis, and this would indeed change the meaning of the physical processes behind the movement in Budyko space.

Budyko-type movement in Budyko space
We define the movement in Budyko space expected from the long-term changes in E0/P (i.e., based only on changes in precipitation and potential evaporation) as "along the Budyko curve".We used the "Budyko-type" analytical climatic model of Yang et al. (2008) expressed by Zhang et al. (2015) as where the parameter n represents the contributing effect of catchment characteristics (e.g., vegetation, soils, topography, seasonality in precipitation and potential evaporation, snow-rain characteristics) in each basin, and where the suffix "b" relates to the Budyko nature of the estimate of E/P by this formulation (E/P)b.We solved the value n for each hydrological basin using the mean values of E/P and E0/P obtained from the CMIP5 models (Table 1) in Equation 4from the periods 1910-1939 and 1961-1990 and calculated (E/P)b for the third and fourth periods, 2010-2039 and 2070-2099.Note that each hydrological basin should have a characteristic Budyko curve as there is one specific value of n for each basin.The direction (θb) and magnitude (Ib) of movement along the Budyko curves is then calculated based on (E/P)b using once again Eq. 2 and 3.
The comparison between movements estimated from CMIP5 projections and those along the Budyko curve can help explain the different drivers behind changes in the evaporative ratio.
For instance, if the CMIP5 movements resemble those along the Budyko curve, we can ratify a significant role of long-term changes in the aridity index as a driver of changes in the evaporative ratio.On the other hand, if CMIP5 movements deviate from the Budyko curves, we expect significant contributions to evaporative ratio change from other drivers such as vegetation and land cover change, seasonality in precipitation and potential evaporation and snow-rain characteristics.We assume no deviation when directions of movement from CMIP5 estimates and along the Budyko curve (θb and θ) fall in the same 10-degree interval.We used Jaramillo and Destouni's (2014) approach to illustrate hydroclimatic change as 'windroses' that summarize change for a large set of river basins.

Results
The selected large river basins cover a wide variety of Earth's hydroclimatic conditions (Fig. 2).The E0/P ranges from river basins where atmospheric energy demand is low and precipitation high─ such as those in Scandinavia or Canada─ to regions experiencing the opposite such in the Sahel and Australia (Fig. 2a).Water partitioning from precipitation into evaporation on the Earth's surface, described by E/P, is also in general higher in the latter river basins than in the former.The distribution of hydroclimatic conditions of the early 20th century  is more skewed towards low E0/P in comparison to future conditions in the late 21 st century (2071-2100), when anthropogenic climate change will have already increased the Earth's temperature by about 4 o C in RCP8.5 (IPCC, 2014) (Fig. 2b).The higher temperatures of the future can preliminarily explain the shift of the distribution of E0/P towards higher aridity (Fig. 2b).During 1911During -1940, 168 , 168 river basins were considered humid since their aridity index was less than two (Barrow, 1992;Greve et al., 2014); by 2071-2100, only 33 river basins will be humid, confirming recent studies finding areas of low E/P to experience a considerable increase in their aridity index (Feng and Fu, 2013;Greve et al., 2019;Lin et al., 2018;Park et al., 2018).1911-1940to 1961-1990 (Fig. 3b (Fig. 3b), striking differences emerge between these change periods and 2011-2040 to 2071-2100 (Fig. 3c).In the latter, movements in Budyko space will converge to predominantly horizontal directions (80 o <θ<100 o ), with a dominant increase in E0/P accompanied by relatively minor changes in E/P across all continents.The few river basins with decreasing E0/P will also move horizontally, such as the Shebelle river basin in Africa and the Krishna, Godavari, Tapti and Mahanadi river basins in India.Thus, it appears that horizontal directions of movement in the space E/P vs E0/P will be the new norm, regardless of the water and energy availability conditions in the river basins and the magnitude of their change.Contrary to the first two change periods, the distribution of directions from 2011-2040 to 2071-2100 is mostly unimodal in terms of E0/P, as most basins experience an increase (Fig. 5e).
In this future change period, as global warming increases temperatures worldwide, the movements along the Budyko curves are almost always horizontal and towards increasing E0/P (70 o <θb<90 o ; Fig. 5f, h), for both RCP4.5 and RCP8.5 scenarios.Nevertheless, the CMIP5 simulations show that only 55% (Fig. 5e) and 63% (Fig. 5g) of river basins will move along their Budyko curves, respectively, leaving 45% and 37% of river basins with movements that deviate from the curves.
For RCP8.5 and the mean of the nine models, the number of river basins where movements from CMIP5 simulations follow the Budyko curves increases from 166 in the first change period to 204 in the third (Fig. 6 and Table 2).The Institut Pierre-Simon Laplace Climate  Lastly, even though CMIP5 projections follow the Budyko curves in 63% of the river basins from 2011-2040 to 2071-2100, there is no significant linear relationship between (E0/P)

Modeling
and (E/P)b (Table 2) as correlation coefficients (R 2 ) are never above 0.04.This applies to all nine Earth system models used.Hence, changes in (E/P)b or movements along the Budyko curve should not be used to predict E/P or movements from CMIP5 simulations.
Figure 6.Deviations of CMIP5-movements (RCP8.5)from their Budyko curves.We assume no deviation when directions of movement from CMIP5 estimates and along the Budyko curve (θb and θ) fall in the same 10-degree interval (See.Fig. 5).

Discussion
We see that movement in the Budyko space of river basins worldwide will converge in the future towards more horizontal directions, regardless of their water and energy availability conditions.The "horizontality" of these movements arises in gains in terrestrial energy availability for evaporation (Arora, 2002;Brutsaert and Parlange, 1998;Donohue et al., 2010) and increasing atmospheric thirst (Falkenmark et al., 2004), accompanied by relatively small changes in the evaporative ratio.As such, the ratio of changes in E/P to changes in E0/P drastically decrease from the second to the third change period, regardless of the model used (Fig. S2).The evaporative ratio will remain largely unaffected despite increasing aridity in basins worldwide.Furthermore, the fact that there is no relationship between changes in the evaporative ratio from CMIP-projections and theoretical Budyko-type estimates (E/P)b (Table 2) calls for caution when predicting changes in runoff or evaporation via the Budyko framework.
The negligible relationship between Δ(E0/P) and Δ(E/P) may also explain why only 41% (1911-1940 to 1961-1990) and 38% (1961-1990 to 2011-2040) of river basins follow their corresponding Budyko curves.Past global studies calculating movement in Budyko space based on precipitation, temperature and runoff observations had found that 74% of a set of almost 900 river basins worldwide were deviating from their Budyko curves in the 20 th century (Jaramillo and Destouni, 2014).Hence, future CMIP5 projections show considerably fewer deviations of hydrological basins from the trajectories of their Budyko curves than historical observations based on runoff, precipitation and temperature.
Furthermore, the number of hydrological basins that are still deviating from their Budyko trajectories under RCP4.5 and RCP8.5 is surprising.The deviations during 1911The deviations during -1940The deviations during to 1961The deviations during -1990The deviations during and 1961The deviations during -1990 to 2011-2040 mainly occur in hydrological basins experiencing increases in E/P that are larger than expected from their Budyko trajectories (i.e., 0 o <θ<60 o ; Fig. 5a, c).On the contrary, during 2011-2040 to 2071-2100, the deviations are due to a slight decrease in E/P even when being subject to an increase in E0/P (i.e., 90 o <θ<100 o ; Fig. 5e, g), which goes against the principles of water and energy availability, which rather expect a mild increase in E/P (i.e., the shape of Budyko curve).Such deviations are evident across all nine models.We now highlight several possibilities explaining the relatively high number of basins deviating from their Budyko curve trajectories.
First of all, land-use change is known to push the long-term movement of river basins in Budyko space beyond the range of slopes given by a typical Budyko-shaped curve (Destouni et al., 2013;Donohue et al., 2007;Renner et al., 2013).Nevertheless, land-use changes would need to explain the large spatial patterns of disagreement beyond these river basins' borders (Fig. 6c).
For instance, all the river basins with the confluence in the North and Baltic Seas evidence a similar disagreement pattern in 2011-2040 to 2071-2100, and a large-scale vegetation conversion is unlikely to explain this pattern alone.
According to Taylor et al. (2012), all CMIP5 models use the same present and future land cover and land-use scenarios based on Hurtt et al. (2011).However, most CMIP5 participating climate models do not provide simulations of forcing due to only land management and vegetation cover change, making it a challenge to ratify if all disagreements of movement in Budyko space are related to these.The Land-Use and Climate, Identification of Robust Impacts (LUCID) experiment found large uncertainties regarding the impacts of vegetation change and land management on hydroclimatic characteristics, based on seven CMIP5 (Brovkin et al., 2013;Noblet-Ducoudré et al., 2012;Pitman et al., 2009) and CMIP6 simualtions (Hurtt et al., 2020).
Our study uses the outputs of three of the models included in the LUCID project; the Max Planck Institute for Meteorology Earth system model (MPI-ESM), the Institut Pierre-Simon Laplace Climate Modeling Center model (IPSL-CM) and the Model for Interdisciplinary Research on Climate (MIROC).Other efforts to disentangle the effects of land use on the hydroclimate based on CMIP5 projections (Kumar et al., 2013) have found that areas with large land cover changes experience a net increase in summer surface albedo, decrease in summer evaporation and increase in summer temperature (i.e., in North America and Eurasia) which may explain some of the disagreements.
Second, adding to the complexity of the effects on land-use change on water partitioning are the potential effects of "greening" of the Earth system by CO2 fertilization (Zeng et al., 2016) or its counterpart; a water-saving response that reduces stomatal conductance and transpiration (Ainsworth and Rogers, 2007;Betts et al., 2007;Medlyn et al., 2001).If any of these two effects could explain a large number of deviations from Budyko-type trajectories, it would be the second, as most deviations exhibit decreases in E/P in a considerable number of basins, despite an increase in E0/P.However, there are also deviations during the first two periods, 1911-1940 to 1961-1990 and 1961-1990 to 2011-2040, when the increase in CO2 emissions is not yet as significant as in the future.Everything depends on how well do CMIP5 models can simulate the delicate balance between the stomatal closure effect of increased atmospheric CO2 on plants and the fertilization effect arising from increasing vegetation or greening (Zeng et al., 2016;Zhang et al., 2016;Zhu et al., 2016;Jaramillo et al., 2018).This improved understanding of vegetation responses under global warming is essential for ecohydrological adaptation and coping strategies under future climate change (Singh et al., 2020).
Thirdly, the deviations of a large number of basins from their Budyko curves may also be attributed to the potential evaporation model chosen to quantify the aridity index.Global studies have seen differences in estimates of aridity when using different potential evaporation models (Greve et al., 2019), such as Penman-Monteith-based models or the radiation-based method here used or its corrections.Nevertheless, it is worth stating that although these imprints some range of uncertainty to the estimates of the aridity index and then of horizontal movement, it is the vertical movement (i.e., Δ(E/P)) that seems to contribute mainly to the deviations of movements from their Budyko-curve trajectories.
Fourthly, long-term intra-annual changes in energy and water availability related to seasonality may also account for deviations from the Budyko curves (Chen et al., 2013;Zanardo et al., 2012); the evaporative ratio may gradually change if precipitation patterns shift within the year, even with the same total annual precipitation.For instance, if precipitation shifts from months of high to low potential evaporation, the amount of precipitation partitioning into actual evaporation will decrease, decreasing the evaporative ratio (Xing et al., 2018).Similarly, a precipitation shift from snow to rain due to higher temperatures in winter and spring will decrease runoff (Berghuijs et al., 2017(Berghuijs et al., , 2014)), which under constant conditions of annual precipitation will increase the evaporative ratio, moving basins upward in Budyko space.
manuscript submitted to Water Resources Research Based on these last points, it appears that the Budyko framework is not enough to represent all hydrological change, which the CMIP5 models may indeed capture.The Budyko framework is based on the spatial distribution of catchments at a specific point in time and not necessarily on the temporal distribution of catchments across time (Budyko, 1974).Furthermore, the non-stationarity of climate parameters and basin characteristics related to water storage, vegetation and land cover pose complications to resolving changes in water fluxes via the framework and separating all drivers of change.
To date, the interpretation of movement in Budyko space has been used for a large set of applications such as to determine: 1) how hydroclimatic change manifests in different biomes (van der Velde et al., 2014), 2) hydroclimatic change global effects of water use and water footprint estimations (Jaramillo and Destouni, 2015;Sun et al., 2021), 3) the influence of forest characteristics on water yield resilience to climate warming (Creed et al., 2014), 4) hydroclimatic change and implications for land water management (Piemontese et al., 2019), 5) the existence of shifts in hydroclimatology (Heidari et al., 2021) and drought (Maurer et al., 2021), and 6) the hydrological effects of vegetation change (Chen et al., 2021).Our study indicates that although we foresee a dominant effect of climate change and global warming in the partitioning of precipitation on land and movement in Budyko space, around half of the largest basins in the world will still deviate from their Budyko curves.Regardless of the reasons for such findings, techniques quantifying and separating the climatic and non-climatic drivers of hydrological change will still be needed for future attribution of changes.In addition, these techniques help quantify human water consumption, water footprint, and impacts of humans on the water cycle and water resources.

Conclusions
We find that 1901 to 2100 movements in Budyko space (in this case, E/P vs E0/P) are predominantly horizontal and will become more horizontal in the foreseeable future.The trend towards horizontal directions of movement arises as global warming increases the aridity index accompanied by much smaller changes in the evaporative ratio.Although the rate of change of movements of most basins increased from 1901 to 2100, this increase results from increasing potential evaporation worldwide and the variables used to construct this space.We find that with global warming, more hydrological basins will also converge towards movement along their Budyko curves; however, 37% will still deviate from their Budyko curve trajectories under RCP8.5.The deviations correspond to a slight decrease in the evaporative ratio and a high increase in the aridity index, which goes against the water and energy availability principles of the Budyko framework and implies less evaporation than expected by the framework.Such deviations can be explained by land-use and vegetation changes or shifts in the seasonality of precipitation or snow to rain ratios.

Figure 1 .
Figure 1.Movement in Budyko space based on CMIP5 projections.The horizontal component

Figure 2 .
Figure 2. The mean values from the nine models for a) aridity index (E0/P; yellow) and evaporation ratio (E/P; purple) during the period 1910 to 2100 for 353 large river basins and during the historic and RCP8.5 scenarios.b) Location in Budyko space of the river basins according to the 30-year periods 1911-1940 (red) and 2071-2100 with RCP8.5 (blue), along with the corresponding Budyko curves based on the mean parameter n of all basins.The water and energy limits constraining water and energy availability are straight black lines.

Figure 3 .
Figure 3. Direction of movement in Budyko space (θ) in the three change periods under the historical and RCP8.5 scenarios.

Figure 4 .
Figure 4. Magnitude of movement in Budyko space (I; Eq. 3) in the three change periods under the historical and RCP8.5 scenarios.

Figure 5 .
Figure 5. Roses of movement in Budyko space for the three change periods based on CMIP5 simulations (a, c, e, g; green roses) and according to the Budyko type model by Yang et al.(2008) (b, d, f, h; grey roses).We construct roses for the three change periods1911-1940 to   1961-1990 (a, b), 1961-1990 to 2011-2040 (c, d)  under the RCP8.5 scenario, and for the change period 2011-2040 to 2071-2100 under both e) RCP4.5 and g) RCP8.5 scenarios.The range of Center model, IPSL-CM5A-MR, and the Norwegian Earth System Model, NorESM1-M, present the highest agreements for the last change period 2011-2040 to 2071-2100.The number of basins with movements following the curves increases into the future as movements become more horizontal in Budyko space, with dominating increases in E0/P and small changes in E/P.Entire regions where basins were deviating from their Budyko curves in the first two change periods will start following the curve trajectories in the future, such as mainland North America (e.g., Mississippi River basin) and the large basins of Mediterranean Europe and Western Asia.Nevertheless, 68 river basins such as the Orinoco, Paraiba do Sul and Tocantins in South America, and Potomac, Saint Lawrence and Savannah River basins in North America consistently deviate from the Budyko curves across the three change periods.South-East Asia is,in particular, a region where river basins consistently deviate.Furthermore, movements in 26 river basins will deviate from their Budyko curves across all nine models, including the Indigirka and Kolima in Asia, the Niger, Orange, Senegal and Gambia in Africa, the Negro and Colorado (Argentina) rivers in South America, or the Back, Coppermine and Yukon Rivers in North America (Suppl.Materials).

Table 1 .
List of CMIP5 climate models analyzed in the study, including their spatial resolutions and main model components.

Table 2 .
Nr. of basins where movement in Budyko space according to CMIP5 projections follows the Budyko curves and statistics of linear reggressions between changes in E/P and (E/P)b for the last change period.