# Budyko Framework Based Analysis of the Effect of Climate Change on Watershed Evaporation Efficiency and Its Impact on Discharge Over Europe

Julie Collignan, Jan Polcher, Sophie Bastin, Pere Quintana-Seguí

► To cite this version:

Julie Collignan, Jan Polcher, Sophie Bastin, Pere Quintana-Seguí. Budyko Framework Based Analysis of the Effect of Climate Change on Watershed Evaporation Efficiency and Its Impact on Discharge Over Europe. Water Resources Research, 2023, 59 (10), pp.e2023WR034509. 10.1029/2023WR034509. insu-03892710v2

HAL Id: insu-03892710

https://insu.hal.science/insu-03892710v2

Submitted on 17 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR034509

**Key Points:**

* Evaluate the role of climate change on the evolution of the discharge of European rivers

* Using the Budyko model as a detection tool for changes in hydrological behaviors of watersheds

* Analyze the effects of changes in the intra-annual distribution of precipitation on the evolution of the annual discharge

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

J. Collignan,
julie.collignan@lmd.ipsl.fr

**Citation:**

Collignan, J., Polcher, J., Bastin, S., & Quintana-Segui, P. (2023). Budyko framework based analysis of the effect of climate change on watershed evaporation efficiency and its impact on discharge over Europe. *Water Resources Research*, 59, e2023WR034509. https://doi.org/10.1029/2023WR034509

Received 5 APR 2023
Accepted 22 AUG 2023

**Author Contributions:**

**Conceptualization:** Julie Collignan, Jan Polcher

**Formal analysis:** Julie Collignan

**Methodology:** Julie Collignan, Jan Polcher, Sophie Bastin

**Resources:** Pere Quintana-Segui

**Supervision:** Jan Polcher, Sophie Bastin

**Writing – original draft:** Julie Collignan

**Writing – review & editing:** Jan Polcher, Sophie Bastin, Pere Quintana-Segui

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Budyko Framework Based Analysis of the Effect of Climate Change on Watershed Evaporation Efficiency and Its Impact on Discharge Over Europe

Julie Collignan¹ ID, Jan Polcher¹ ID, Sophie Bastin², and Pere Quintana-Segui³

¹Laboratoire de Météorologie Dynamique/IPSL, Ecole Polytechnique/CNRS, Paris, France, ²Laboratoire Atmosphères, Observations Spatiales/IPSL, CNRS, Paris, France, ³Observatori de l'Ebre, Universitat Ramon Llull, CSIC, Roquetes, Spain

**Abstract** In the context of climate change, the stakes surrounding water availability are rapidly intensifying. Decomposing and quantifying the effects of climate on discharge allows us to understand their impact on water resources better. We propose a methodology to separate the effect of change in the annual mean of climate variables from the effect of the intra-annual distribution of precipitation. It combines the Budyko framework with land surface model (LSM) outputs. The LSM is used to reproduce the behavior of 2,134 reconstructed watersheds across Europe between 1902 and 2010, with climate inputs as the only source of change. We fit a one-parameter approximation of the Budyko framework to the LSM outputs. It accounts for the evolution of the annual mean in precipitation (P) and potential evapotranspiration (PET). We introduce a varying parameter in the equation, representing the effect of long-term variations in the intra-annual distribution of P and PET. To better assess the effects of changes in annual means or intra-annual distribution of P, we construct synthetic forcings fixing one or the other. European results show that the trends in the annual averages of P dominate the trends in discharge due to climate. The second main climate driver is PET, except over the Mediterranean area, where changes in intra-annual variations of P have a higher impact on discharge than trends in PET. Therefore, the effects of changes in the intra-annual distribution of climate variables are to be addressed when looking at changes in annual discharge.

**Plain Language Summary** Water availability is a challenge for all of society. Various competing activities rely on this resource, and its scarcity can lead to social, economic, and environmental conflicts. With climate change, river discharge and, more generally, the full water cycle is impacted. Furthermore, multiple human actions such as dams and irrigation concurrently change the balance of the water cycle over watersheds. To comprehensively understand the dynamics of discharge, it is essential to analyze the potential influence of direct human activities alongside the impacts of natural climatic factors. Models are a way to represent reality with an understanding of the physical phenomena included. They can be used to represent the behavior of watersheds without human intervention. In light of this, we have developed a methodology to highlight the climate factors impacting discharge. Annual discharge changes are driven mainly by changes in annual precipitation over Europe. The increasing temperature leads to an ever-growing evaporative demand and is the second most impacting factor over most of Europe. However, in the Mediterranean area, where water is more limited, changes in the seasonality of precipitation have a higher impact than changes in the evaporative demand.

## 1. Introduction

Water is a key resource for the whole of society and both its excess and scarcity can lead to challenging economic, environmental, and social issues. Understanding the hydrological cycle and how it evolves due to a changing climate is a significant challenge of this century.

Over the past century, several studies have shown the impact of climate change on climate variables in Europe. Annual precipitation increased between 1901 and 2005 over most of Europe except the Mediterranean area, where they tended to decrease (Christidis & Stott, 2022; Douville et al., 2021; Knutson & Zeng, 2018). Trends per decade are less significant due to the high inter-annual variability of precipitation *P* (Douville et al., 2021). Trends in potential evapotranspiration (PET) are linked to an increase in the energy available at the surface, which is highly correlated to rising temperatures (Douville et al., 2021; Vicente-Serrano et al., 2014). Few studies have directly examined European PET trends, except over the Mediterranean area, where studies have shown a
significant increase in PET over the end of the century (Kitsara et al., 2013; Vicente-Serrano et al., 2014, 2019). The intra-annual variations of climatic variables are more difficult to assess, and only a few indices exist to measure the inter-annual changes in the distribution of climate variables. For example, García-Barrón et al. (2018) defined indices to assess the evolution of the intra-annual cycle of *P* over time throughout the Iberian Peninsula. At the end of the century, they identified a shift of the main rainfall periods toward autumn, especially over the Atlantic basins, and an increase in the inter-annual variability of the intra-annual cycle, especially over the Mediterranean basins. For precipitation, studies have shown that not only is the annual average of *P* changing, but there are differences between summer and winter, depending on the area (Christidis & Stott, 2022; Zveryaev, 2004). Moreover, over the past few decades, extreme precipitation events affecting this area have significantly intensified (Ribes et al., 2019). Therefore, it is important to investigate the effects of changes in the annual averages of climate variables along with the effect of changes in seasonality and intra-annual distribution of these variables. The distribution of *P* within the year and its coupling or decoupling from the atmospheric demand PET will influence water partitioning between evapotranspiration and discharge on the annual scale.

Transformations in different climate variables governing the water cycle alter the equilibrium in the water balance over the different watersheds, thus impacting the discharge of rivers. Milly et al. (2005) showed that worldwide discharge trends are and will continue to be significantly impacted by changes in climatic factors. Over Europe, statistically significant trends in discharge are observed in historical records (positive in the northern region and negative in the south and east). These trends are spatially coherent with precipitation changes (Stahl et al., 2010; Vicente-Serrano et al., 2019). Y. Yang et al. (2018) show that discharge is less sensitive to PET changes than to changes in *P*.

The effects of intra-annual variations of *P* on discharge are primarily considered in the literature through the study of seasonality and annual extremes of *P* and PET in order to examine their impact on floods (Douville et al., 2021; Milly et al., 2002; Rottler et al., 2020), drought events (Douville et al., 2021; Vicente-Serrano et al., 2014), and more generally, on discharge peaks (Bouwer et al., 2008; Tuel et al., 2022). Stahl et al. (2010) found the trends in discharge over the end of the century were disconnected between summer flows and winter flows for an ensemble of small near-natural catchments in Europe. Blöschl et al. (2019) showed that increasing autumn and winter rainfall led to increased floods in northwestern Europe.

However, rivers are also highly managed, and human activities are an important driver of change in how watersheds function (Ficklin et al., 2018; Riedel & Weber, 2020). A significant difficulty in analyzing the effect of climate on historic discharge changes is decomposing the effects of the different drivers of change and isolating them from each other to better understand their relative importance (Stahl et al., 2010; Vicente-Serrano et al., 2019). Several studies have concentrated on catchments that are regarded as near-natural or unimpaired in order to investigate the effects of climatic changes on discharge (Stahl et al., 2010; Y. Yang et al., 2018). However, this highly limits the areas studied (Vicente-Serrano et al., 2019), especially in Europe, where the high population density and long history of water management limit the study to small catchments (Stahl et al., 2010).

Another approach is to use models to separate the factors involved in discharge changes; different types of models have been developed. Models relying on a few variables and adjusted parameters are favored for their simplicity of use and interpretation. One example is using statistical models fitted over specific areas, such as linear regressions (Bouwer et al., 2008; Ficklin et al., 2018; Vicente-Serrano et al., 2019). More complex models integrate nonlinear relationships and physical boundaries. However, all these parsimonious models are empirical: they rely on adjusted parameters over the area and the time period studied and lack a clear physical meaning. The parameters often cannot be generalized and transposed to other areas or future climates (Coron et al., 2014; Reaver et al., 2020).

Other methods are physical-based hydrological and land surface models (LSMs). They require more data and computational power. They do not always accurately represent a whole real hydrological system depending on which processes are included in them but allow a meaningful assessment of hydrologic aridity (Y. Yang et al., 2018). They have grown increasingly complex and are able to reproduce the behavior of watersheds and to model “natural flow” regimes (Decharme et al., 2019; Gudmundsson et al., 2017; Schneider et al., 2013; F. Wang et al., 2018). However, due to their complexity, it is more difficult to decompose the effects of individual climate factors and to interpret their outputs than with other simpler models.

In light of this, we propose a tool that combines the simplicity of the more empirical model with the heightened performance and complexity of the physical-based model to better understand the phenomena encapsulated behind the adjusted parameters.
We use here the well-known and widely used empirical Budyko framework (Mianabadi et al., 2020). It is predicated upon utilizing the annual mean of water and energy balances at the watershed scale (Tian et al., 2018), taking into account the water and energy limitations of the physical system. It was initially conceived over multiple catchments. Parametric equations were developed to introduce an empirical parameter adjusting the framework to the specific evaporation efficiency of each catchment over an equilibrium period (H. Yang et al., 2008; L. Zhang et al., 2004). However, equilibrium disruptions, due to climate change or any other direct human activities and vegetation change, highlight limitations to the model. Moreover, most disruptive features are concurrent. The parameter introduced has no evident physical meaning and is just a well-adjusted proxy to E/P over a specific catchment and period. There is no straightforward method to attribute changes in the adjusted parameters to specific climatic or nonclimatic features (Berghuijs et al., 2020; Reaver et al., 2022), as for any parsimonious model with calibrated parameters.

To focus on the effects of climate change, the present study applies the Budyko framework to the outputs of a state-of-the-art LSM. The latter represents the constant physical behavior of watersheds. The only source of change in the dynamics of the modeled watersheds is the evolving climate variables introduced. Using LSM outputs also allows for adjusting the near-surface atmospheric variables to more adequately decompose the effects of the different elements of climate variability and change.

This article is organized as follows: Section 2 covers the methodology developed. It describes the Budyko framework with its underlying hypothesis and limitations and the state-of-art LSM. Then, we describe how we apply the framework to the chosen LSM. Here, we create Synthetic forcings to test if our methodology yields an optimal analysis of the effects of different aspects of climate change. We also explain the use of the time-moving window to examine temporal trends in the different climatic effects. In Section 3, we present the results of the effect of different elements of climate change across Europe (changes in annual averages against changes in the intra-annual distribution of climate variables) on discharge trends over the past century. Section 4 provides a comprehensive analysis of the advantages of our findings, while also highlighting the present constraints and areas for further investigation. Finally, in Section 5, we summarize our conclusions.

## 2. Methodology

### 2.1. The Budyko Framework

#### 2.1.1. General Presentation

Over watersheds considered as closed systems, the water balance Equation 1 applies when explaining the equilibrium between the variables of the hydrological cycle: the river discharge (Q), the evapotranspiration (E), the precipitation (P) and the change in the water storage over the watershed between two-time steps (ΔS).

$$P - \Delta S = Q + E \qquad (1)$$

Long-term, ΔS can be negligible. Ideally, this hypothesis should be applied over a long enough period that the system's equilibrium is reached (L. Zhang et al., 2008). It also supposes no external disturbances impact the water budget, such as groundwater mining or water transfers to or from other basins.

The Budyko framework, which is frequently used in hydrological research to study the partitioning of P into E and Q, draws from this long-term equilibrium of water balance over a catchment coupled with the energy balance. It postulates that the partition of the annual water budget between runoff and evapotranspiration over catchments, represented by the evapotranspiration E, is a function of the relative water supply (rainfall P) and the atmospheric water demand (PET) (Tian et al., 2018; Xing et al., 2018; D. Yang et al., 2007). The latter depends on both available energy and aerodynamic resistance (Barella-Ortiz et al., 2013). Therefore, this framework considers the system's water and energy limitations, which cannot evaporate more than the atmospheric demand allows and more water than the catchment receives from the water source (P). In short, it defines the “Budyko space” (Berghuijs et al., 2020; Reaver et al., 2022).

This framework relies on a closed water budget in time and space, neglecting ΔS. Therefore, it must be applied over a closed watershed and fitted on a long-term equilibrium. To be freed from seasonal water storage variations, we use a time series of a yearly resolution (hydrological year) in this study. For the region considered, the hydrological year starts in September, at the end of the dry season, when the reservoirs are supposedly at their lowest. It
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Budyko framework: relationship between evapotranspiration ratio (E/P) and aridity index (PET/P) (Fu's equation). E, PET, P are annual averages. ω associated with the purple curve is larger than ω associated with the orange curve and translates into a higher evaporation efficiency above the watershed. For a given watershed with constant characteristics, there is still a dispersion around the curve of the dots for a given year due to intra-annual variations of the climate cycle (orange dots). The curve and its associated ω represent the average behavior of the watershed. The framework includes trends in annual climate variables by a displacement along the curve (red arrow). However, it does not include trends that could impact the way water is partitioned over the catchment such as long-lasting trends in the intra-annual distribution of P and PET (blue arrows).</figcaption></figure>

minimizes the differences in ΔS from year to year. Later on, unless specified otherwise, the variables P, E, and Q represent the annual averages over the hydrological year. We then apply the framework over minimum periods of 11 years, considered a long enough period for ΔS to be negligible over most catchments, dependent on the area (Han et al., 2020). We tested this hypothesis with the outputs of the LSM, and we found that ΔS is about a hundred times smaller than Q when 11-year sub-periods are considered (not shown).

### 2.1.2. One Parameter Equation

The original Budyko framework was empirically constructed over a set of catchments to define a curve followed, on average, by catchments in the Budyko space. Different analytical approximations to this hypothesis (Budyko curves) have been developed, expressing the evapotranspiration rate (E/P) as a function of the aridity index (PET/P) over a catchment (Figure 1).

More specifically, the framework was extended to analyze individual catchments over a stable period. Parametric equations were developed which introduced an empirical parameter representing the specific position of the catchment within the Budyko space (H. Yang et al., 2008).

Two of the most widely used are the Fu equation (Equation 2) (Ning et al., 2019; Simons et al., 2020; L. Zhang et al., 2004, 2008; Zheng et al., 2018) and the Mezentsev-Choudhury-Yang equation (Equation 3) (Luo et al., 2020; Roderick & Farquhar, 2011; W. Wang et al., 2020; Xing et al., 2018; Xiong et al., 2020; H. Yang et al., 2008). These can be found under different names in the literature, such as the Tixeront-Fu equation for Equation 2 or Turc-Mezentsev for Equation 3 (Andréassian & Sari, 2019).

$$ \frac{E}{P} = 1 + \frac{PET}{P} - \left(1 + \left(\frac{PET}{P}\right)^{\omega}\right)^{\frac{1}{\omega}} \quad (2) $$

$$ E = \frac{P * PET}{(P^n + PET^n)^{\frac{1}{n}}} \quad (3) $$

The two parameters derived from Equations 2 and 3 are linearly correlated, implying that both equations are almost equivalent (Andréassian & Sari, 2019; Du et al., 2016; Roderick & Farquhar, 2011; H. Yang et al., 2008). We examine the sensitivity of the results to the parametric equation used. We obtain very similar results for the methodology with either equation used. We conclude that we could use either equation. For the rest of the study, we use results obtained with Fu's equation (Equation 2).

E measurements are not available over large spatial and temporal scales. Therefore, most studies work from the analysis of Q, which can be calculated from the water balance Equation 1, where ΔS has been neglected. With Fu's equation (Equation 2) used to express E in Equation 1, it yields (Equation 4):

$$ Q = P * \left(1 + \left(\frac{PET}{P}\right)^{\omega}\right)^{\frac{1}{\omega}} - PET = f(P, PET, \omega) \quad (4) $$

### 2.1.3. Discussion of the Watershed Parameter

The watershed parameter is empirical; it is obtained by fitting data from a specific catchment during a period of assumed equilibrium state. It determines the position of the catchment in the Budyko space.

The specificity of the parameter relates to all factors impacting the evaporation efficiency of the watershed other than changes in the average aridity index (Donohue et al., 2012; Padrón et al., 2017; L. Zhang et al., 2004). The most common hypothesis is that it reflects the various hydrological characteristics of the watershed, such as topography, vegetation coverage, and soil properties, which play a part in the annual partitioning of P into E and Q over the catchment (Gudmundsson et al., 2017; Reaver et al., 2022). Some are considered time-invariant
(soil type, topography, etc.), while others are possibly affected by long-lasting changes. These can occur in the hydrological properties of the surface water system, most likely due to direct anthropogenic activities such as river management, irrigation, and land cover changes. It leads to the "catchment trajectory conjecture" (Reaver et al., 2022), which suggests that the watersheds would follow an average Budyko-curve (Figure 1, red arrow) if it were not for changes in hydrological properties independent of changes in the average aridity index.

Several studies attempted to analyze the evolution of watershed behavior between two equilibrium states, a period of reference and a period of post-changes (Jiang et al., 2015; Luo et al., 2020; W. Wang et al., 2020; Zhao et al., 2018; Zheng et al., 2018) and then fit the parameter independently over each period. Two distinct curves (Figure 1) were acquired using distinct watershed parameters to characterize the pre- and post-change equilibrium states. As a first hypothesis, they then considered that deviation from the initial curve (period of reference) is only due to changes in the land surface, such as the effect of anthropogenic activities and land cover variations. Assuming $\omega$ to be climate invariant, the changes due to climate are considered in the framework only through the modifications of the average P/PET (Figure 1, red arrow). It follows the hypothesis that watersheds follow their Budyko curve if the catchment's surface characteristics remain constant.

However, studies have shown that not all catchments under climate change exhibit this behavior. There is a climate dependence of the deviation to the initial curve. Reaver et al. (2022) showed that reference catchments with the long-term stability of land use did not always follow their Budyko-curve. With the previous hypothesis, this deviation could be misinterpreted as a change in the land-surface characteristics. Padrón et al. (2017) found that the variability in the parameter is highly correlated to climate features such as snow fraction precipitation and the storm arrival rate. Over their extensive global database, the correlation between vegetation indices and direct anthropogenic influence factors is only secondary. Jaramillo et al. (2022) used CMIP6 multi-models ensemble to fit Budyko curves over several basins for the period 1901–1950 and to calculate ET/P PET/P for 2051–2100. They compared the results of the ensemble to those obtained with the hypothesis that catchments should follow their initial Budyko curve. Most basins will not follow the curve under climate change, showing a climate dependence of the deviation from the initial curve.

To circumvent the limitation due to the hypothesis of $\omega$ being climate invariant, several studies have tried to locate an expression of the watershed parameter as a function of pertinent factors. It would allow us to express the evolution of $\omega$ over time and decompose the effects of climate and human activities through the different factors chosen. If valid, it would also allow transposing the expression to unmonitored catchments where $\omega$ cannot be directly fitted or to future catchment conditions. Different methods, such as step-wise regressions and neural networks, were used to identify pertinent factors. Such methods require enough information on the chosen factors; strong hypotheses stand behind the expression.

Most studies construct their function across several basins, accounting for spatially different human, climate, and land characteristics (D. Li et al., 2013; S. Li et al., 2022; Ning et al., 2019; Tian et al., 2018; Xing et al., 2018; X. Zhang et al., 2019). A variety of factors were selected: environmental factors such as soil moisture, seasonality of P and PET, aridity index (S. Li et al., 2022; Ning et al., 2019), vegetation fraction and routing depth (Gentine et al., 2012; D. Li et al., 2013; Ning et al., 2019), relief ratio, drought severity index, seasonality of P, and synchronicity between P and PET (Xing et al., 2018); direct human factors such as irrigated surfaces (Tian et al., 2018), the amount of water applied for irrigation (D. Li et al., 2013), land use, land cover change in highly managed areas (Tian et al., 2018), and even gross domestic product per capita (X. Zhang et al., 2019). The chosen factors are highly dependent on the area studied.

Another strong hypothesis is that such a relationship defined over spatial differences is applicable to explain temporal differences (Luo et al., 2020). Other studies (Jiang et al., 2015; Zhao et al., 2018) looked at time-varying human activities and climate change to construct expressions, using a time-moving window to fit the evolution of the catchment parameter over a basin. This approach faces another limitation due to the availability of information on the different factors' time evolution. Ning et al. (2019) used a mixed technique, applying their fit across 30 basins at different time scales using moving time windows and found that the impact of vegetation cover and climate seasonality on the watershed parameter was stronger over longer time steps, showing that the weight of different factors varies with the time scale. Ultimately, the pertinent factors highly differ among studies, regions, time periods, and climate types (Padrón et al., 2017).

Moreover, recent studies question the hypothesis underlying these studies and "the catchment trajectory conjecture" (Berghuijs et al., 2020; Reaver et al., 2020, 2022). The study demonstrates that the parameter exhibits a lack of independence from climate but also depends on the biophysical characteristics of the catchment directly due to
the dependence of *E*, PET, or *P* on those features (Reaver et al., 2020). The highly nonlinear relationships between all the features involved in the evaporation efficiency of the catchment and the average *P* and PET contradict the hypothesis that the parameter *ω* can be expressed as a simple function of independent parameters. It also explains why previous studies were so different from one another. The catchment parameter is, therefore, a mathematical tool to represent the evaporation efficiency of a catchment over a given period and has no physical meaning in itself. It is not transferable through time and space. It only positions the catchment within the Budyko space (Reaver et al., 2022). It can still be used to study the position of the catchment in the Budyko space. Studying the deviation from the curve may provide insight into how factors besides aridity affect the water balance (Berghuijs et al., 2020).

## 2.2. Simulations With a LSM

To isolate the climate change effect from other factors that could affect watersheds, we work with the outputs of a LSM. The model constructs watersheds with constant hydrological properties and represents an idealized watershed without any direct changes from human activities and other nonclimatic disturbances. Therefore, the only source of long-term change would be due to a difference in response to an evolving climate.

### 2.2.1. A “Natural Reference” Simulation

This study uses the LSM Organizing Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) from the Institut Pierre Simon Laplace. It includes biophysical and biogeochemical processes to simulate the global carbon cycle and quantify terrestrial water and energy balance. It runs coupled to an atmospheric model or in stand-alone conditions with an independent data set to force the atmospheric conditions. Here, we use the model in stand-alone conditions, forced with the data set GSWP3 covering 1901–2013 (Hyungjun, 2017) at the resolution of 0.5° for all climate variables.

The hydrological network of the ORCHIDEE LSM is constructed from the hydrological elevation model HYDROSHEDS (Lehner et al., 2008), which covers the area studied with the resolution of 30 arc seconds (approximately 1 km at the equator). The hydrological information is then upscaled to the resolution of the atmospheric grid, the hydrological coherence being preserved by the construction of hydrological transfer units at the subgrid level (Polcher et al., 2022). From a database of gauging stations, upstream basins are reconstituted on the hydrological elevation model grid and then projected on the atmospheric grid during the process. We have access to 2,134 stations over the area studied for which the LSM calculates a discharge and for which we have the reconstituted upstream basin (Figure S1 in Supporting Information S1).

The LSM ORCHIDEE, more specifically the Schématisation des EChanges Hydriques à l'Interface Biosphère-Atmosphère module, uses the USEB (unstressed surface energy balance) method to model PET (details in Barella-Ortiz et al., 2013). This method relies on the turbulent diffusion equation to calculate the potential soil evaporation PET<sub>soil</sub>, obtained from the air density, the aerodynamic resistance, and the humidity gradient. The USEB method estimates the virtual surface temperature from an unstressed surface-energy balance, computing a new energy balance considering a saturated surface (Barella-Ortiz et al., 2013). Potential transpiration is driven by PET<sub>soil</sub> but limited by vegetation resistance, calculated in LSM ORCHIDEE and based on plant functioning types maps and LAI calculations (Guimberteau et al., 2012). Then, the total PET is calculated by summing the potential evaporation and the potential transpiration. PET is reduced to the actual evapotranspiration *E* by a “moisture availability function” (Barella-Ortiz et al., 2013).

Over the course of several years, the model has been tested and validated on many aspects of the land surface processes (hydrology, vegetation, and carbon cycle processes). This attests to its quality to reproduce the water and energy balance and also discharge over different areas over the globe (Guion et al., 2022; Nguyen-Quang et al., 2018; Polcher et al., 2022; Tafasca et al., 2020; F. Wang et al., 2018). Comparing the LSM outputs directly to observations for discharge is challenging, mainly due to the absence of certain processes in the models, including those resulting from direct human activities and the extensive water and river management (F. Wang et al., 2018), as it is the case in our area of interest. Based on previous literature, we can assume the model proficiently emulates the mechanisms underlying actual evaporation, thereby effectively replicating the “natural” response of watersheds with persistent physical attributes to the past climate conditions prevalent in Europe. We study *Q* variations and not the absolute value of *Q* since we know that the output of the LSM does not represent the complete processes over real catchments. We focus here on the impact of the changes in atmospheric parameters
<table><caption>**Table 1**
Synthetic Forcings Created</caption><thead><tr><th></th><th>Forcing name</th><th>Average <i>P</i></th><th>Intra-annual variation of <i>P</i></th><th>Description<sup>a</sup></th></tr></thead><tbody><tr><td>1</td><td><i>ref</i></td><td>—</td><td>—</td><td>Reference forcing: GSWP3 (1901–2012)</td></tr><tr><td>2</td><td><i>f</i> 2000</td><td>Fixed</td><td>Fixed</td><td><i>P</i> has been entirely fixed for each year, equal to the precipitation and the seasonality of the year 2000.</td></tr><tr><td>3</td><td><i>cstmean</i></td><td>Fixed</td><td>—</td><td>Only the average value of <i>P</i> has been fixed for every year to the one of the year 2000</td></tr><tr><td>4</td><td><i>cstintravar</i></td><td>—</td><td>Fixed</td><td>Only the intra-annual variations of <i>P</i> have been fixed for every year to the one of the year 2000</td></tr></tbody></table>

on land surface responses with constant characteristics. The modeled watersheds react to the climate data input at each time step (30 min time step). Therefore, the LSM output depends on both the evolving annual average and the evolving distribution over the year of the climate variables.

For consistency in the calculation of *E* and PET, we take both from the output of ORCHIDEE forced with GSWP3. The gridded outputs (PET, *E*) are at the resolution of the forcing data set (0.5°). *P* is the sum of rainfall and snowfall in GSWP3. Then, we consider the annual mean *P*, PET, and *E* over hydrological years, integrated over each catchment. The catchments' shape has been reproduced at a finer resolution and then projected on the 0.5°grid.

### 2.2.2. Synthetic Forcings to Analyze the Effect of Variation of Seasonality

In order to better understand the effect of inter and intra-annual climate variations on the Budyko framework and on discharge *Q*, we construct synthetic climate forcings, fixing one or the other.

The calculation of PET includes many related climate variables and nonlinear relationships, making it very difficult to anticipate how a change in a given climate variable may influence its behavior. It is, therefore, too complicated to create synthetic forcings for which we can modify climate variables to fix PET seasonality, for instance. Therefore, we only modify the precipitation *P* in the synthetic forcings to see how it impacts our results compared to the reference forcing.

The reference forcing is the GSWP3 data set from September 1901 to September 2012 (3 hr time step). Then we constructed three forcings, which were modified over hydrological years (Table 1, Figure 3a):

* *f* 2000: A forcing where all 3 hr values of *P* are set to the values of the year 2000 (September 1999 to September 2000) for each year. Therefore, all components of *P* (average and intra-annual variations) are set constant.

* *cstmean*: A forcing for which we keep the relative intra-annual distribution of *P* of each year, but where the average *P* of each year is set constant. The 3 hr values of *P* are scaled so the hydrological year average is set to the one of the year 2000 (September 1999 to September 2000).

* *cstintravar*: A forcing for which we keep the annual average of *P* for each year, but where the relative intra-annual distribution of *P* is set constant. The 3 hr values of *P* are set to the values of the year 2000 (September 1999 to September 2000) for each year and then scaled over each hydrological year so the yearly average is set to the one of the corresponding years in the reference forcing.

### 2.3. Combining Both Models

In this study, we apply the Budyko framework to the output of an LSM to explore the sensitivity of the empirical parameter to climate change and the resulting effect on discharge. The watersheds in the LSM have constant

<footnote><sup>a</sup>For forcings 2–4, *P* has been modified compared to the reference: the average value of *P* over the year and/or the distribution of precipitation over the year have been fixed for each year to the value of the year 2000.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Scheme of the method: the land surface model (LSM) is obligated with the forcing data set to calculate *E*. The LSM is considered to represent the “climatic reality” over a catchment without any changes in the watershed characteristics. We then average *P*, PET, and *E* and integrate them over each watershed to get annual averages for all catchments. Then, we fit the Fu equation. (a) The fit of the equation over the entire century results in the calculation of an empirical parameter $\bar{\omega}$, which represents the average catchment characteristics. (b) To have an evolution of $\omega_t$ over time, the fit was then successively applied over an 11-year sliding time period.</figcaption></figure>

biophysical characteristics. The LSM then reproduces the interaction of the land surface with climate parameters. It is affected by climate change and no other source of change. We use the changes in $\omega$ as a proxy for changes in the partitioning of *P* into *Q* and *E* other than direct changes in average PET/P; it focuses on the deviation from the initial curve and attempts to decompose its dependence on climate. In this case, any deviation to the curve is only due to climate effects. Since $\omega$ has no clear physical meaning, we do not analyze directly the changes in $\omega$ but rather how they impact the evolution of discharge.

Using an LSM, we can also change various climate parameters to better address how they weigh in the modeled changes. We develop a varying $\omega_t$ to capture part of the change in the evaporation efficiency of the watersheds due to climate. We compare its effects to the magnitude of change in discharge already captured with the traditional framework, which only considers changes in annual averages of PET/P.

### 2.3.1. Fit of the Evaporation Efficiency Parameter $\omega$

The watershed parameter of the Budyko curve is calculated over each catchment with a fit of the equation curve $E/P = f(PET/P)$ (Equation 2), using the minimum root mean square error for a given set of annual averages of evapotranspiration *E*, precipitation *P* and PET data (Jiang et al., 2015; D. Yang et al., 2007). We fit the parameter once with all points over the entire period covered by the climate data set to obtain $\bar{\omega}$ representing the average behavior for each catchment (Figure 2a).

For a watershed with constant hydrological properties (which is the case when considering modeled watershed in ORCHIDEE), if we consider the “catchment trajectory conjecture,” $\omega$ is independent of climate, and the catchment follows its initial curve. However, $\omega$ varies for a given watershed because of climate. For instance, over an equilibrium state, intra-annual variations of the climate cycle induce a variability of the annual values (E/P, PET/P) around the fitted curve. The distribution of rain changes the covariance between *P* and PET over the year. A difference in storm depth over a catchment can change the capacity of the soil to store water, the response of vegetation, and change the dynamic of the water partition into runoff and evaporation even if the annual amount
(a) Example for give watershed. Map: the colors show the share of the grid point within the watershed. Yellow points are completely within it while bordering grid points are red. Modified forcings over a given basin: the first row of graphs shows the inter-annual variability of *P* (cyan), *PET* (green) and *E* (purple) for each forcings. The second row shows the average seasonal distribution of *P* over the catchment for each forcings over the entire century. The average monthly distributions of *P* shown here are computed over the century. It however varies from one year to another for the reference forcing and the forcing *cstmean* which is not illustrated here.

6217110: Guadalquivir, Rio : Cantillana, Lon: -5.83° Lat: 37.59°

<figure><img src="image_5.png" /></figure>

(b) Watershed parameter $\bar{\omega}$ fitted over the entire time period (dashed purple line) and $\omega_t$ fitted successively over a sliding 11-year time-window (blue line) for the reference forcing.

<figure><img src="image_6.png" /></figure>

(c) Discharge estimated with Budyko for the reference forcing: $Q = f(P, PET, \omega_t)$ (blue line), $Q_c = f(P, PET, \bar{\omega})$ (purple line), $Q_\omega = f(P_{rand}, PET_{rand}, \omega_t)$ (black line) with their associated trends. Unsignificant trends are dashed. Here all trends are significant.

<figure><img src="image_7.png" /><figcaption>Figure 3.</figcaption></figure>
of precipitation stays constant (Donohue et al., 2012). More generally, a change in synchronization between P (water available) and PET (energy demand from the atmosphere) will change E/P for the same average climate (Abatzoglou & Ficklin, 2017; S. Li et al., 2022). In an equilibrium state, the intra-annual variations should be without trends and only result in white noise around that equilibrium. The fitted parameter ω represents the average behavior of the basin. For a catchment under climate change, however, this variability could lead to a significant permanent deviation from the initial curve if this intra-annual distribution tends to have a trend.

### 2.3.2. Introducing a Varying Watershed Parameter ωᵢ

With its simple framework, the Budyko model does not cover possible changes at intra-annual time scales. The average effect of this synchronization is included in the adjustment parameter ω, which is, therefore, not completely independent of climate. Therefore, long-term changes in seasonality should induce a climatic time dependence which is not accounted for in the framework with a constant ω. Therefore, considering a varying parameter should improve the Budyko model to reproduce E/P and its climatic evolution.

To obtain a varying parameter ωᵢ for each catchment, we carry out several fits over successive 11-year time-sliding sub-periods (Figure 2b). We chose 11 years as the smallest time length to apply the Budyko framework relevantly, considering that each 11-year sub-period is stationary (ΔS = 0). This allows us to focus on long-term changes and to minimize the impact of year-to-year “transient” effects (e.g., soil storage and groundwater changes) (Y. Yang et al., 2018). Tian et al. (2018) found that below a certain time length, the fit of the ω parameter was too unstable to be relevant.

### 2.3.3. Decomposing the Impact of Climate on Discharge Trends

The watershed parameter ω is a conceptual variable that provides little insight into the magnitude of discharge changes. Thus, we examine the impact of ωᵢ on the river discharge Q and compare these changes to the impact of annual averages of climate variables (P and PET) changes on Q over time. To simplify the discussion, we gather the annual averages of P and PET in a “climate” variable C = (P, PET).

Following our previous hypothesis (Equation 4), Q can be estimated with the Budyko framework using C and ω: $Q = f(C, ω)$.

Q can be decomposed with first-order partial derivatives (Equation 5), with the first term of the right-hand side representing the partial derivative due to climate variables C and the second term for the partial derivative due to changes in the watershed parameter ωᵢ. We then estimate the partial derivatives due to C and due to ω independently.

$$ \frac{dQ}{dt} = \frac{\delta Q}{\delta C} \frac{dC}{dt} + \frac{\delta Q}{\delta \omega} \frac{d\omega}{dt} \quad \text{with } C = (P, PET) \quad (5) $$

To independently estimate the partial derivative due to climate variables C, we must cancel the second term (Equation 5, left side). To do so, we calculate the discharge $Q_c = f(C, \bar{\omega})$, with a constant value of ω. The trend of that discharge $\frac{dQ_c}{dt}$ matches the term with the partial derivative due to C in Equation 5.

To estimate the partial discharge trend due to ωᵢ, we need to eliminate the trends in annual averages of P and PET over the century to cancel the first term (Equation 5, left side). We randomly draw P and PET pairings for each year. We do so several times and average the results for each year. It gives us a random climate without trends over the century. We then apply Fu's equation (Equation 2) with the resulting random annual averages of P and PET and the varying ωᵢ calculated with the forcing before the random drawing. It gives $Q_ω = f(C_{rand}, ωᵢ)$ for which the climate trends are only due to variations captured by the time-varying parameter ωᵢ. The trend $\frac{dQ_ω}{dt}$ matches the term with the partial derivative due to ω in Equation 5. In the end, we get:

**Figure 3.** Time series obtained through the full application of our methodology for a given basin in Spain. (a) shows the inter-annual variability of annual averages of climate variables P, potential evapotranspiration (PET), and E modeled by the land surface model, forced with the different synthetic forcings. E mostly relates to P. (b and c) are results for the reference forcing. (b) shows the varying ωᵢ resulting from the time-sliding window calculation (blue curve), compared to $\bar{\omega}$ calculated with one fit over the entire century (dashed purple line). (c) shows the decomposition of the discharge, comparing the full discharge to partial discharges and their respective trends. The full discharge Q is modeled with Fu's equation with annual averages of P and PET from the reference forcing and ωᵢ. The first partial discharge $Q_C$ is the one calculated with the constant parameter $\bar{\omega}$. It covers most Q variations for the given basins. The second partial discharge $Q_ω$ covers some of the missing variations of Q and some of the missing trends due to deviations to the average curve. From that figure, we can conclude that most variations and trends of the discharge in this basin are explained by C = (P, PET).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Boxplot of Nash-Sutcliffe coefficient (NSC) for all watersheds: for the forcing of reference with the constant parameter $\bar{\omega}$, with the varying parameter $\omega_t$ and for the forcing cstintravar (where the seasonal distributions of P have been fixed over the entire time period) with a constant $\bar{\omega}$. It represents how well the Budyko model reproduces the discharge output from Organizing Carbon and Hydrology In Dynamic Ecosystems. A value above 0.5 is considered satisfactory. Very similar results are found when looking at $R^2$ from a linear regression.</figcaption></figure>

$$ \frac{dQ}{dt} = \frac{dQ_c}{dt} + \frac{dQ_{\omega}}{dt} \qquad (6) $$

We calculate the trends of each term and their significance using the Mann-Kendall nonparametric test, associated with the Thiel-Sen slope estimator. It gives us time series and associated trends for each studied watershed. Figure 3 shows an example of a watershed in southern Spain.

# 3. Results

## 3.1. Performance of Budyko With or Without a Variant Parameter $\omega$

We hypothesize that for watersheds with constant hydrological properties, the dispersion of annual points around the curve is due to intra-annual variations of climate. If these variations did not exist, catchments would follow their Budyko curve, and we could use it to model the discharge almost perfectly.

To test this hypothesis, we examine the performance of the Budyko curve with a constant parameter $\bar{\omega}$ to reproduce the discharge from the LSM for the reference forcing compared to the forcing cstintravar. For that latter forcing, we removed the intra-annual variations of P from 1 year to another, which should render the performance of the Budyko curve model close to perfect if the hypothesis is valid.

We use the Nash-Sutcliffe coefficient (NSC) as a performance indicator (Equation 7, Figure 4). We consider an NSC >0.5 to be satisfactory (Moriasi et al., 2007).

$$ NSC = 1 - \frac{\sum_{i=0}^{\text{years}} (Ql_i - Qb_i)^2}{\sum_{i=0}^{\text{years}} (Ql_i - \bar{Q})^2} \quad \begin{cases} \text{with } Ql = \text{dicharge from the LSM} \\ \text{and } Qb = \text{Result from the methodology with Fu's equation} \end{cases} \qquad (7) $$

We obtain NSC values above 0.5 for 89.9% of all 2,134 watersheds tested for the Budyko curve with a constant parameter ($Q_c$, calculated with a constant $\bar{\omega}$) applied with the reference forcing (boxplot on the left, Figure 4). Therefore, the average curve model is rather effective in reproducing the annual discharge over watersheds with constant hydrological properties reacting to an evolving climate.
For the forcing *cstintravar*, NSC for $Q_c$ increases to above 0.6 for all watersheds (boxplot on the right, Figure 4). It confirms our hypothesis: the average Budyko curve model is even more effective if there are no intra-annual variations of *P* from 1 year to another. Therefore, most of the variability that is not captured by the average Budyko curve over the past century is due to the intra-annual variability of *P* and the covariance of *P* and PET.

When looking at NSC for the framework applied to the reference forcing with a varying parameter $Q(\text{ref}) = f(C(\text{ref}), \omega_t)$, we gain up to 0.26 points of NSC for the tested watershed and reach 94.1% of all watersheds with NSC >0.5 (boxplot on the center, Figure 4). It does not reach the performance to reproduce $Q_c$ with the forcing *cstintravar*. However, it allows to catch some of the deviation to the curve due to intra-annual trends of climate variables. We capture long-term trends following our choice of the 11-year time-moving window. It validates our hypothesis that introducing a varying watershed parameter $\omega_t$ improves the framework to better encompass climate variability and the effect of climatic trends on discharge, including the effect of climate change on the intra-annual distribution and covariance of climate variables (*P* and PET).

To sum up, for watersheds with constant hydrological properties under historical climate, most of the deviation to the average curve model (i.e., changes in the evaporation efficiency of catchments) is due to variations in the intra-annual distribution of climate variables (*P* and PET). Our varying parameter improves the framework by allowing us to capture the long-term trends of these variations. We now analyze their effect on the discharge and compare them to the direct effect of trends in the annual average of climate variables.

## 3.2. Comparing the Effects of Intra-Annual Variations of *P* on Discharge *Q* to the Effects of Variations in Annual Averages of *P* in Europe

We consider our area of study, western Europe (2,134 watersheds modeled) (Figure 5). To better illustrate our results, we also take two contrasted basins: one in Italy (Figure 6) and another in England (Figure 7).

Figures 5a–5c show the relative trends over each basin for the reference forcing, respectively, of $Q$, $Q_c$, and $Q_\omega$. There are significant decreases in the total discharge $Q$ (Figure 5a) (–0.3% to –0.4% per year over the past century) over sparse basins in Spain, the Pyrenees, Italy, Slovenia, Greece, and Eastern Europe. There are significant increases (Figure 5a) (+0.2% to +0.4% per year over the past century) over sparse basins in France, Germany, Denmark, Sweden, Northern UK, and Serbia. These trends are primarily due to changes in the annual averages $C = (P, \text{PET})$ since the average Budyko curve model $Q_c$ captures most of the signal (Figure 5b). The inter-annual variability of $C$ is high, making the trends less than 95% significant over most basins for $Q$ and $Q_c$. Both selected catchments better illustrate it (Figures 6b and 7b): for the reference forcing (top left), the dominant effect in the variations of annual discharge $Q$ (blue line) is due to the annual mean of climate variables $C$ (purple line). Clearly, both curves have very similar high inter-annual variations and trends.

Changes in $C$ are the dominant factors explaining the climatic trends in $Q$ over the past century in Europe. The results obtained with the forcing *cstintravar* (bottom right for Figures 6b and 7b and maps Figures 5j–5l) confirm it. It shows that without inter-annual changes in *P* distribution (in other words, with a maximum reduction of the inter-annual changes in the annual covariance of *P* and PET), the discharge *Q* obtained and the associated relative trends are very similar to the results obtained with the reference forcing. Therefore, the effects of changes in the annual covariance of *P* and PET are minor compared to the effects of changes in the annual mean of climate variables *C* in most of Europe.

However, in some areas, the effects of the intra-annual distribution of *P* should be addressed. If we look at the Tiber River in Italy (Figure 6b), the trend in $Q_c$ (purple line) is significant for both the reference forcing and the forcing *cstintravar*. However, the total discharge $Q$ (blue line) trend is only significant for the forcing *cstintravar*. For the reference forcing, the decreasing trend in the discharge due to $C(Q_\omega)$ is counteracted by the increasing trend due to changes in the evaporation efficiency ($Q_\omega$), making the final trend in discharge $Q$ insignificant.

More generally, over Europe, when we erase the inter-annual variability of $C$, we capture the effect of trends in the intra-annual distribution of *P* and PET, through changes in the evaporation efficiency, in $Q_\omega$ (Figure 5c). It tends to increase discharge, especially in southwestern Spain, Italy, and the west of France (+0.1% per year over the century). It corresponds to the increasing trend of the black line in Figure 6b, top left graphs for the Tiber River. It has an opposite trend toward decreasing discharge in eastern Europe and has a relatively neutral effect in the rest of the continent (Figure 5c and, in the example of the English basin, Figure 7b, top left graph, black line). It amplifies the trends due to annual averages $C$ changes over certain watersheds such as the Duero basin (north-western Spain, decrease in discharge), western France, and northern Germany. Indeed, we note a significant increase in discharge over certain watersheds where the effect of changes in $C$ alone was insignificant.
(a) Relative trends $\frac{dQ/\bar{Q}}{dt}$: ref

<figure><img src="image_9.png" /></figure>

(b) Relative trends $\frac{dQ_c/\bar{Q}}{dt}$: ref

<figure><img src="image_10.png" /></figure>

(c) Relative trends $\frac{dQ_\omega/\bar{Q}}{dt}$: ref

<figure><img src="image_11.png" /></figure>

(d) Relative trends $\frac{dQ/\bar{Q}}{dt}$: f2000

<figure><img src="image_12.png" /></figure>

(e) Relative trends $\frac{dQ_c/\bar{Q}}{dt}$: f2000

<figure><img src="image_13.png" /></figure>

(f) Relative trends $\frac{dQ_\omega/\bar{Q}}{dt}$: f2000

<figure><img src="image_14.png" /></figure>

(g) Relative trends $\frac{dQ/\bar{Q}}{dt}$: cstmean

<figure><img src="image_15.png" /></figure>

(h) Relative trends $\frac{dQ_c/\bar{Q}}{dt}$: cstmean

<figure><img src="image_16.png" /></figure>

(i) Relative trends $\frac{dQ_\omega/\bar{Q}}{dt}$: cstmean

<figure><img src="image_17.png" /></figure>

(j) Relative trends $\frac{dQ/\bar{Q}}{dt}$: cstintravar

<figure><img src="image_18.png" /></figure>

(k) Relative trends $\frac{dQ_c/\bar{Q}}{dt}$: cstintravar

<figure><img src="image_19.png" /></figure>

(l) Relative trends $\frac{dQ_\omega/\bar{Q}}{dt}$: cstintravar

<figure><img src="image_20.png" /></figure>

(m) Colorbar: significant relative trends (% of change per year over the century)

<figure><img src="image_21.png" /><figcaption>Figure 5.</figcaption></figure>
In other areas, such as the Tiber River in Italy, or in southern UK, the intra-annual variability of *P* and PET counteracts the effect of *C*, making the relative total *Q* trends lose their significance due to opposite signals. We note the decreasing trend is due to the evolution of *C*, while the effect of the change in the intra-annual distribution of the climate variables tends to increase the discharge.

In order to investigate the impacts of intra-annual variations of *P* on discharge, we analyze the results of the synthetic forcing *f* 2000 and *cstmean* (respectively top right and bottom left Figures 6b and 7b and maps Figures 5d–5f and 5g–5i). For the synthetic forcing *f* 2000 (Figures 5d–5f), *P* have been entirely set for each year to *P* of the year 2000. Therefore, this only yields the trends due to changes in PET, both for changes due to annual climate variables and changes in the evaporation efficiency of the catchment. For the synthetic forcing *cstmean*, only the annual mean of *P* has been set. In this case, the trends are due to PET and changes in the intra-annual distribution of *P*.

For the forcing *f* 2000, the effect of PET is toward a decrease in discharge over all of Europe (less than −0.1% to −0.2% per year over the century) (Figure 5d). For both the chosen examples, the effect of PET (top right graphs) tends to decrease discharge (purple line, *Q*<sub>c</sub> when *P* has been fixed). It is consistent with the significant increase in PET (Figures 6a and 7a, top right). The effect of intra-annual variations of PET on changes in the evaporation efficiency (Figure 5f and black lines, top right graph Figures 6b and 7b) has the same order of magnitude, if not a little smaller (less than −0.1% per year over the century), than the effect of inter-annual change of the annual average of PET (Figure 5e or purple line top right graph Figures 6b and 7b). It tends to amplify the latter's effect, especially over western France and southern UK. It has a slightly opposite effect toward increasing trends in *Q* (less than +0.08% per year over the century) over the east of Europe, west of Spain, and for the Tiber river. The effect of changes in the annual mean of PET, in this specific case, is canceled in the total discharge (blue line) by the effect of the changes in the intra-annual distribution of PET captured in *Q*<sub>ω</sub> (black line) (Figure 6b).

For the forcing *cstmean*, we now add the effect of changes in the intra-annual covariance of *P* and PET due to changes in the intra-annual distribution of *P*. Depending on the area, there are two different responses. The two basins chosen in the example each correspond to one type of response. In the case of the basin in England (Trent River), the results obtained for the forcing *cstmean* (Figure 7b, bottom left) are very similar to the results obtained for *f* 2000 (Figure 7b, top right). This means that the effect is due to changes in intra-annual synchronicity of *P* and PET has little impact compared to the effect of the annual mean of PET over that particular basin. It matches the results over northern Europe, especially over France, Germany, and southern UK, where the trends in *Q* (Figure 5g) are mainly driven by changes in the annual mean of PET (Figure 5h). However, over the Tiber River in Italy, the results obtained for the forcing *cstmean* (Figure 6b, bottom left) shows that the changes in the total discharge *Q* (blue line) match the changes due to the evolution of *ω*<sub>t</sub> (*Q*<sub>ω</sub>, black line). In this latter case, the effect of the intra-annual variations of *P* is dominant compared to the effect of changes in PET. This matches the results over southern Europe (Spain, Italy) where for the forcing *cstmean*, the trends in *Q* (Figure 5g) are driven mainly by changes in the evaporation efficiency (Figure 5i). This increase in discharge diverges from the trends due to changes in *C* in the area (reference forcing and forcing *f* 2000, purple lines).

The discharge trends for both forcings, namely *f* 2000 and *cstmean*, are statistically significant across multiple watersheds, independent of the high inter-annual variability observed in the annual mean of *P*. Trends are significant for 1,883 basins with the forcing *f* 2000 and 1,756 for the forcing *cstmean* against only 352 basins with significant trends in *Q* out of 2,134 for the reference forcing. However, the magnitude of these trends is also quite small. Comparing the discharge obtained with the reference forcing shows that the main factor driving *Q* is the annual mean of *P* since the discharge trends look entirely different, when free of its variations.

To sum up, the results obtained with the synthetic forcings, the annual mean of *P* is the first driver of change in the annual discharge over all of Europe. However, its high inter-annual variability tends to hide the trends in most areas. The second most important climatic driver of discharge change depends on the area. Over southern Europe (Italy, Spain), where water is the limiting factor to evapotranspiration, the second most important climatic factor driving discharge changes is the intra-annual distribution of *P*. Over the rest of Europe, where water is less limiting, the second most important factor driving discharge changes is the increasing PET.

**Figure 5.** Decomposition of significant relative *Q* trends (% of change per year over the century) for all the tested forcings: the first line is the reference forcing. The first column is the total change in *Q*, the second is the partial change due to trends in the annual average of *P* and potential evapotranspiration (PET), and the last column is the partial change due to changes in the watershed parameter, mostly due to trends in the intra-annual distribution of *P* and PET. For the modified forcings: *f* 2000 has the annual average and intra-annual distribution of *P* fixed for every year to their value for the year 2000. *cstmean* has only the annual average of *P* fixed. *cstintravar* has only the intra-annual distribution of *P* fixed. White areas do not have significant trends.
Tiber River:
Roma

<figure><img src="image_22.png" /><figcaption>(a) Annual average of climate variables *P* (light blue line), *PET* (yellow line) and *E* (purple line) modeled with the LSM for each academic forcing. Not shown here, the intra-annual distribution of *P* has been fixed for the forcings *f2000* and *cstintravar*.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Discharge estimated with Budyko and their respective trends: full discharge modeled $Q = f(P, PET, \omega_t)$ (blue line), discharge $Q_c = f(C, \bar{\omega})$ (purple line) with only variations in $C$ accounted for, and $Q_\omega = f(P_{rand}, PET_{rand}, \omega_t)$ (black line) with only the variations of evaporation efficiency, mostly due to the intra-annual covariance of $P$ and $PET$, accounted for. Unsignificant trends are dashed. Here the results for each academic forcing are shown. We represented the normalized anomaly of discharge $((Q - Q_{mean})/Q_{mean})$ in order to better compare the plots to each other. The scale of the y-axis changes and is divided by 5 for the forcings $f2000$ and $cstintravar$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 6.** Example 1: time series obtained through the full application of our methodology for a given basin in Italy.</figcaption></figure>
Trent River:
Colwick

<figure><img src="image_25.png" /></figure>

(a) Annual average of climate variables *P* (light blue line), *PET* (yellow line) and *E* (purple line) modeled with the LSM for each academic forcing. Not shown here, the intra-annual distribution of *P* has been fixed for the forcings *f2000* and *cstintravar*.

<figure><img src="image_26.png" /></figure>

(b) Discharge estimated with Budyko and their respective trends: full discharge modeled $Q = f(P, PET, \omega_t)$ (blue line), discharge $Q_c = f(C, \bar{\omega})$ (purple line) with only variations in $C$ accounted for, and $Q_\omega = f(P_{rand}, PET_{rand}, \omega_t)$ (black line) with only the variations of evaporation efficiency, mostly due to the intra-annual covariance of $P$ and $PET$, accounted for. Unsignificant trends are dashed. Here the results for each academic forcing are shown. We represented the normalized anomaly of discharge $((Q - Q_{mean})/Q_{mean})$ in order to better compare the plots to each other. The scale of the y-axis changes and is divided by 5 for the forcings *f2000* and *cstmean*.

<figure><img src="image_27.png" /><figcaption>**Figure 7.** Example 2: time series obtained through the full application of our methodology for a given basin in England.</figcaption></figure>
4. Discussion

Similar to the results of several studies (Abatzoglou & Ficklin, 2017; Jaramillo et al., 2022; S. Li et al., 2022; Padrón et al., 2017; Reaver et al., 2022; Xing et al., 2018), we find that the average Budyko curve model with a constant watershed parameter $\bar{\omega}$ does not capture climate-related changes in the watershed behavior impacting its evaporation efficiency. Even with constant hydrological land surface characteristics, most catchments do not follow their average curve over the past century. The deviation to the curve can significantly affect *Q*'s long-term trends over the past century if we free our analysis from the high inter-annual variability of *P*. It is in accordance with the results of Reaver et al. (2022) criticizing the “catchment trajectory conjecture.”

Parameter $\omega$ has no direct physical meaning but is a proxy to represent the evaporation efficiency of catchments (Berghuijs et al., 2020; Reaver et al., 2022). However, since it cannot be expressed as a function of clearly defined factors, it is difficult to attribute the changes in the evaporation efficiency to specific climatic features (Berghuijs et al., 2020). Y. Yang et al. (2018) assume that further reductions in *Q* declining trends due to changes in catchment properties are likely associated with elevated atmospheric CO₂ concentration or increased rainfall intensity. Other studies find correlations between changes in the evaporation efficiency of catchments and storm depth, the portion of precipitation such as snow (Donohue et al., 2012; Padrón et al., 2017; Xing et al., 2018). Using the outputs of an LSM, our studies allow us to test a selection of hypotheses by adjusting climate parameters. We find that the climatic deviation to the average Budyko curve over the historical records is mainly due to variations in the intra-annual distribution of *P*.

We introduce a time-varying window to fit the parameter of Fu's equation in order to capture trends in the deviation to the average curve in the Budyko space. The choice of window size determines the temporal scale of accounted-for trends. This functions as a frequency filter and only captures the effect of variations over periods the size of the window or larger. We must balance the length of our data set and the appropriate length of the trends we choose to analyze. Since our aim is to investigate the effects of climate change, we do not need to capture the high inter-annual variability and can focus on decadal trends or longer (Y. Yang et al., 2018). Furthermore, a shorter time window would not be adapted to the hypothesis of the Budyko framework, which needs a long enough period to be fitted. So, the window cannot be shorter (Tian et al., 2018). An exploration of an extended time window could be conducted to investigate the limited duration of time that captures the most significant impact on discharge. However, the longer the time window, the fewer points we will have to evaluate the trends.

In our methodology, we decompose the trends due to climatic changes in evaporation efficiency and the trends due to changes in average climate variables *P* and PET. One limitation in our decomposition method is that the variations in the evaporation efficiency captured in the deviation to the average curve are not entirely independent from the variations of average *P* and PET. The relationship between *P*, PET, and evaporation efficiency is complex and relies on many interrelated factors (Reaver et al., 2020). We find in our study that the changes captured in the varying $\omega_t$ are mostly due to changes in the covariance of the intra-annual distribution of *P* and PET. However, the effect of the intra-annual distribution of climate variables on discharge is not completely independent from the annual mean of *P* and PET because of the difference in sensitivity of the system to a change in water availability. It can impact the magnitude of the identified trends. It is shown by the slight differences observed in $Q_\omega$ between the reference forcing and the forcing cstmean (Figures 5c and 5i) and between the forcing f 2000 and the forcing cstintravar (Figures 5f and 5l). For each pair of forcings, the intra-annual distribution of *P* is the same, but the inter-annual mean of *P* differs. The difference in $Q_\omega$ for each pairing is due to a link between the annual mean and the intra-annual distribution of *P*. Therefore, the amplitude of the effect of the intra-annual distribution of *P* and PET quantified here may depend on the choice of the fixed average *P* (again, *P* from the year 2000 in this study). The observed differences were found to be comparatively insignificant in light of the identified trends, indicating that the fundamental findings regarding Europe would remain unchanged; therefore, we opted to disregard them. When studying specific basins, it could be interesting to choose specific pairings of intra-annual distributions/ annual averages of *P* to construct synthetic forcings, to compare how specific associations combine.

Furthermore, we cannot simply fix PET or its intra-annual variations in our synthetic forcings due to its nonlinearity dependence on a number of climate variables. Therefore, we are unable to decompose the effects of PET as easily as for the effects of *P*, which would be interesting to do, especially in the areas where *P* is less limiting, such as in western France or northern Europe.
Our methodology allows for the separation of the effect of primary and secondary climatic drivers on discharge trends. We look at the trends in *P*, and PET for the forcing GSWP3. Our results concur with those in the literature, validating that this forcing reasonably reproduces the climatic trends of the past century over Europe. The trends in PET are significantly (95% level) increasing over Europe. However, the trends in *P* are most often nonsignificant because of its high inter-annual variability, with a significant trend in the annual average of *P* for 413 catchments out of 2,134 selected. The present study finds that the main driver of annual discharge *Q* (trends and inter-annual variability) is the annual mean of *P*. As expected with the increase in *P* over western Europe and the decrease in *P* observed in the Mediterranean area (Christidis & Stott, 2022; Douville et al., 2021; Knutson & Zeng, 2018), the trends in *Q* have followed the same direction. It concurs with the finding of Stahl et al. (2010) and Vicente-Serrano et al. (2019), who found strong spatial consistency between streamflow changes and global rainfall changes.

H. Yang et al. (2008) show that *Q* is universally more sensitive to changes in *P* than to changes in PET, for a fixed land surface condition. Similarly, we find that over most of Europe, the second most important climatic factor on discharge changes is PET, which leads to a decrease in discharge due to the increasing evaporative demand by the atmosphere. Over the Iberian Peninsula and the Mediterranean area, however, PET trends have a lesser impact. There, the water limit is the prevailing factor, having been attained by the end of spring and persisting throughout the entirety of summer. Therefore, a warmer summer does not have a strong impact. The evolution of intra-annual variations of *P* is the second most important factor impacting the changes in the annual discharge, with a higher effect on discharge than the increase of PET over the past century. The intra-annual covariance of *P* and PET impacts the annual behavior of the catchment and the annual balance between evapotranspiration and discharge since it changes the timing between water and energy available throughout the year. The evolution of the intra-annual cycle of *P* tends towards decreasing discharge in the Mediterranean area. It partially counteracts the effect of decreasing *P* and increasing PET on discharge. Therefore, the intra-annual distribution of *P* deserves more attention when studying the evolution of annual discharge. In most studies, it is only considered to look at changes in discharge peaks, floods, or droughts (Douville et al., 2021; Milly et al., 2002; Rottler et al., 2020; Tuel et al., 2022; Vicente-Serrano et al., 2014). We calculate the indices defined by García-Barrón et al. (2013) to evaluate the trends in the intra-annual cycle of *P* for the forcing GSWP3. Similarly to the authors' findings, in Spain, we identify a shift over the end of the century towards a more bimodal distribution of precipitation throughout the year. However, the trends in the intra-annual cycle are mostly qualitative. The tendencies of the annual cycle to have an increasingly marked seasonality, concentrating rain events in fewer but more extreme events over the year, can explain the increasing runoff and relative discharge. Our methodology allows to identify these effects despite the only qualitative trends observed in the indices that measure the intra-annual distribution of *P*.

We apply our parametric model to LSM outputs to isolate the discharge variations due to changes in climate factors. This methodology relies on the capacity of the chosen LSM to reproduce the “natural” response of a catchment to climate, such as its behavior and response to changes in the intra-annual distribution of *P*. The amplitude of our results could depend on the choice of the LSM or the forcing data. We tested the use of other forcing data sets: WFDEI (Weedon et al., 2014), which covers the period from 1979 to 2010, with the same resolution as GSWP3, and E2OFD (Beck et al., 2017), while also covering 1979–2010 but at a lower resolution. We also tested another model, SURFEX (Quintana-Seguí et al., 2020), forced with SAFRAN (Quintana-Seguí et al., 2017), over the Ebro river. This yielded similar results over the overlapping period with little differences in the trends' significance and amplitude. This indicates that the resolution of the forcing exerts a greater influence on the results compared to a specific forcing or model employed. This confirms the suitability of utilizing an LSM as a climatic reference in accordance with our methodology. In the future, when looking at specific basins, it would be interesting to use higher resolution forcings to obtain a more accurate picture of the effects of climate change on discharge. In this case, the diversity of behaviors exhibited among subbasins within a given catchment could be elaborated upon by distinguishing the behavior of upstream subbasins within mountainous regions from that of the downstream portion, which may display differential responses to climate change.

## 5. Conclusion

Our methodology combines a physical-based model to a parsimonious model. The first allows to identify the climatic changes in the empirical parameter of the second. The second allows for a simple decomposition of the relative changes in discharge. In this case, the Budyko framework and a one-parameter equation: the deviation
from the average curve corresponds to a change in the evaporation efficiency of the catchment. The state-of-the-art LSM was used to simulate changes in the evaporation efficiency under the climate of the past century, independent from any other disruptive process. The successive fit of the parametric equation allows us to find the climatic dependence of the deviation to the average curve in the Budyko space over time.

For a given catchment, we quantify its effect on *Q*. Over the past century, the primary climatic source of deviation to the average curve is the change in the intra-annual distribution of *P*. We compare the impact of that deviation on changes in the average annual discharge compared to the change due to average climate variables *P* and PET. Over Europe for the past century, the main climatic driver of change in the average *Q* is the change in the average *P*. The second main driver of discharge change is PET over most of Europe except the Mediterranean area, where a change in the intra-annual distribution of *P* weighs more on *Q* changes than PET. Therefore, the effect of the intra-annual distribution of *P* should be addressed when studying the evolution of the average discharge and water availability under climate change, especially over the Mediterranean.

If we were to work from observations instead of model outputs, there would be other non-climate-related sources of variability, such as direct human activities or vegetation changes which would modify watershed behavior. Our next step is to apply the methodology to quantify these human-induced changes and compare their magnitude to those attributed to climate change in the present study's responses.

## Data Availability Statement

The forcing data set GSWP3 used to grid *P* and other climate data and run the LSM over Europe between 1901 and 2010 in the study is freely available upon registration (Hyungjun, 2017). The LSM used to calculate PET and model the discharge in this study is ORCHIDEE (IPSL [Institut Pierre Simon Laplace], 2017), available on their website. The outputs used for this study at the annual time step for each catchment are gathered in a file freely available on Zenodo.org (Collignan et al., 2023). Stations used in the study come from the Global Runoff Data Centre (GRDC) (2022), completed with the Geoportal of Spain Ministerio (Ministerio para la Transición Ecológica y el Reto Demográfico, 2020) and over France with data from the database HYDRO (Ministère de l'écologie, du développement durable et de l'énergie, 2021), where the data are freely accessible but have to be gathered region by region and station by station. The file on Zenodo.org (Collignan et al., 2023) also includes the list of the stations used in the study and their main related metadata: their location and the size of the upstream area used to position the station on the grid. The upstream watersheds are reconstructed using the hydrological elevation model HydroSHEDS (Lehner et al., 2008) to construct the routing graphs for rivers on the LSM grid.

## Acknowledgments

We would like to acknowledge the support of the Agence Nationale de la Recherche under contract HLIAISE (ANR-19-CE01-0017-02). The lead author would like to thank Institut Polytechnique de Paris for the Gaspard Monge fellowship which funded her Ph.D. thesis.

## References

* Abatzoglou, J. T., & Ficklin, D. L. (2017). Climatic and physiographic controls of spatial variability in surface water balance over the contiguous United States using the Budyko relationship. *Water Resources Research*, 53(9), 7630–7643. https://doi.org/10.1002/2017WR020843

* Andréassian, V., & Sari, T. (2019). Technical note: On the puzzling similarity of two water balance formulas—Turc–Mezentsev vs. Tixeront–Fu. *Hydrology and Earth System Sciences*, 23(5), 2339–2350. https://doi.org/10.5194/hess-23-2339-2019

* Barella-Ortiz, A., Polcher, J., Tuzet, A., & Laval, K. (2013). Potential evaporation estimation through an unstressed surface-energy balance and its sensitivity to climate change. *Hydrology and Earth System Sciences*, 17(11), 4625–4639. https://doi.org/10.5194/hess-17-4625-2013

* Beck, H. E., van Dijk, A. I. J. M., Levizzani, V., Schellekens, J., Miralles, D. G., Martens, B., & de Roo, A. (2017). MSWEP: 3-hourly 0.25° global gridded precipitation (1979–2015) by merging gauge, satellite, and reanalysis data. *Hydrology and Earth System Sciences*, 21(1), 589–615. https://doi.org/10.5194/hess-21-589-2017

* Berghuijs, W. R., Gnann, S. J., & Woods, R. A. (2020). Unanswered questions on the Budyko framework. *Hydrological Processes*, 34(26), 5699–5703. https://doi.org/10.1002/hyp.13958

* Blöschl, G., Hall, J., Viglione, A., Perdigão, R. A. P., Parajka, J., Merz, B., et al. (2019). Changing climate both increases and decreases European river floods. *Nature*, 573(7772), 108–111. https://doi.org/10.1038/s41586-019-1495-6

* Bouwer, L. M., Vermaat, J. E., & Aerts, J. C. J. H. (2008). Regional sensitivities of mean and peak river discharge to climate variability in Europe. *Journal of Geophysical Research: Atmospheres*, 113(D19), D19103. https://doi.org/10.1029/2008JD010301

* Christidis, N., & Stott, P. A. (2022). Human influence on seasonal precipitation in Europe. *Journal of Climate*, 35(15), 5215–5231. https://doi.org/10.1175/JCLI-D-21-0637.1

* Collignan, J., Polcher, J., Bastin, S., & Quintana-Segui, P. (2023). Output of the land surface model ORCHIDEE over river catchments in Europe, run with GSWP3 and synthetic forcings where the precipitation is modified [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.8211025

* Coron, L., Andréassian, V., Perrin, C., Bourqui, M., & Hendrickx, F. (2014). On the lack of robustness of hydrologic models regarding water balance simulation: A diagnostic approach applied to three models of increasing complexity on 20 mountainous catchments. *Hydrology and Earth System Sciences*, 18(2), 727–746. https://doi.org/10.5194/hess-18-727-2014

* Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J.-P., Alias, A., et al. (2019). Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications. *Journal of Advances in Modeling Earth Systems*, 11(5), 1207–1252. https://doi.org/10.1029/2018MS001545
Donohue, R. J., Roderick, M. L., & McVicar, T. R. (2012). Roots, storms and soil pores: Incorporating key ecohydrological processes into Budyko's hydrological model. *Journal of Hydrology*, **436**(437), 35–50. https://doi.org/10.1016/j.jhydrol.2012.02.033

Douville, H., Raghavan, K., Renwick, J., Allan, R., Arias, P., Barlow, M., et al. (2021). Water cycle changes. In *Climate change 2021: The physical science basis* (pp. 1055–1210). Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. https://doi.org/10.1017/9781009157896.010

Du, C., Sun, F., Yu, J., Liu, X., & Yaning, C. (2016). New interpretation of the role of water balance in an extended Budyko hypothesis in arid regions. *Hydrology and Earth System Sciences*, **20**(1), 393–409. https://doi.org/10.5194/hess-20-393-2016

Ficklin, D. L., Abatzoglou, J. T., Robeson, S. M., Null, S. E., & Knouft, J. H. (2018). Natural and managed watersheds show similar responses to recent climate change. *Proceedings of the National Academy of Sciences*, **115**(34), 8553–8557. https://doi.org/10.1073/pnas.1801026115

García-Barrón, L., Aguilar-Alba, M., Morales, J., & Sousa, A. (2018). Intra-annual rainfall variability in the Spanish hydrographic basins. *International Journal of Climatology*, **38**(5), 2215–2229. https://doi.org/10.1002/joc.5328

García-Barrón, L., Morales, J., & Sousa, A. (2013). Characterisation of the intra-annual rainfall and its evolution (1837–2010) in the southwest of the Iberian Peninsula. *Theoretical and Applied Climatology*, **114**(3), 445–457. https://doi.org/10.1007/s00704-013-0855-7

Gentine, P., D'Odorico, P., Lintner, B., Sivandran, G., & Salvucci, G. (2012). Interdependence of climate, soil, and vegetation as constrained by the Budyko curve. *Geophysical Research Letters*, **39**(19), 19404. https://doi.org/10.1029/2012GL053492

Global Runoff Data Centre (GRDC). (2022). In situ river discharge [Dataset]. Retrieved from https://www.bafg.de/GRDC/EN/02_srvcs/21_tmsrs/riverdischarge_node.html

Gudmundsson, L., Greve, P., & Seneviratne, S. I. (2017). Correspondence: Flawed assumptions compromise water yield assessment. *Nature Communications*, **8**(1), 14795. https://doi.org/10.1038/ncomms14795

Gudmundsson, L., Seneviratne, S. I., & Zhang, X. (2017). Anthropogenic climate change detected in European renewable freshwater resources. *Nature Climate Change*, **7**(11), 813–816. https://doi.org/10.1038/nclimate3416

Guimberteau, M., Laval, K., Perrier, A., & Polcher, J. (2012). Global effect of irrigation and its impact on the onset of the Indian summer monsoon. *Climate Dynamics*, **39**(6), 1329–1348. https://doi.org/10.1007/s00382-011-1252-5

Guion, A., Turquety, S., Polcher, J., Pennel, R., Bastin, S., & Arsouze, T. (2022). Droughts and heatwaves in the Western Mediterranean: Impact on vegetation and wildfires using the coupled WRF-ORCHIDEE regional model (RegIPSL). *Climate Dynamics*, **58**(9), 2881–2903. https://doi.org/10.1007/s00382-021-05938-y

Han, J., Yang, Y., Roderick, M. L., McVicar, T. R., Yang, D., Zhang, S., & Beck, H. E. (2020). Assessing the steady-state assumption in water balance calculation across global catchments. *Water Resources Research*, **56**(7), e2020WR027392. https://doi.org/10.1029/2020WR027392

Hyungjun, K. (2017). Global soil wetness project phase 3 atmospheric boundary conditions (Experiment 1) [Dataset]. Data Integration and Analysis System (DIAS), 5. https://doi.org/10.20783/DIAS.501

IPSL (Institut Pierre Simon Laplace). (2017). Organizing Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) [Software]. Retrieved from https://orchidee.ipsl.fr/

Jaramillo, F., Piemontese, L., Berghujs, W. R., Wang-Erlandsson, L., Greve, P., & Wang, Z. (2022). Fewer basins will follow their Budyko curves under global warming and fossil-fueled development. *Water Resources Research*, **58**(8), e2021WR031825. https://doi.org/10.1029/2021WR031825

Jiang, C., Xiong, L., Wang, D., Liu, P., Guo, S., & Xu, C.-Y. (2015). Separating the impacts of climate change and human activities on runoff using the Budyko-type equations with time-varying parameters. *Journal of Hydrology*, **522**, 326–338. https://doi.org/10.1016/j.jhydrol.2014.12.060

Kitsara, G., Papaioannou, G., Papathanasiou, A., & Retalis, A. (2013). Dimming/brightening in Athens: Trends in sunshine duration, cloud cover and reference evapotranspiration. *Water Resources Management*, **27**(6), 1623–1633. https://doi.org/10.1007/s11269-012-0229-4

Knutson, T. R., & Zeng, F. (2018). Model assessment of observed precipitation trends over land regions: Detectable human influences and possible low bias in model trends. *Journal of Climate*, **31**(12), 4617–4637. https://doi.org/10.1175/JCLI-D-17-0672.1

Lehner, B., Verdin, K., & Jarvis, A. (2008). New global hydrology derived from spaceborne elevation data [Dataset]. Eos, Transactions American Geophysical Union, 89(10), 93–94. https://doi.org/10.1029/2008EO100001

Li, D., Pan, M., Cong, Z., Zhang, L., & Wood, E. (2013). Vegetation control on water and energy balance within the Budyko framework. *Water Resources Research*, **49**(2), 969–976. https://doi.org/10.1002/wrcr.20107

Li, S., Du, T., & Gippel, C. J. (2022). A modified Fu (1981) equation with a time-varying parameter that improves estimates of inter-annual variability in catchment water balance. *Water Resources Management*, **36**(5), 1645–1659. https://doi.org/10.1007/s11269-021-03057-1

Luo, Y., Yang, Y., Yang, D., & Zhang, S. (2020). Quantifying the impact of vegetation changes on global terrestrial runoff using the Budyko framework. *Journal of Hydrology*, **590**, 125389. https://doi.org/10.1016/j.jhydrol.2020.125389

Mianabadi, A., Davary, K., Pourreza-Bilondi, M., & Coenders-Gerrits, A. M. J. (2020). Budyko framework; towards non-steady state conditions. *Journal of Hydrology*, **588**, 125089. https://doi.org/10.1016/j.jhydrol.2020.125089

Milly, P. C. D., Dunne, K. A., & Vecchia, A. V. (2005). Global pattern of trends in streamflow and water availability in a changing climate. *Nature*, **438**(7066), 347–350. https://doi.org/10.1038/nature04312

Milly, P. C. D., Wetherald, R. T., Dunne, K. A., & Delworth, T. L. (2002). Increasing risk of great floods in a changing climate. *Nature*, **415**(6871), 514–517. https://doi.org/10.1038/415514a

Ministere de l'ecologie, du developpement durable et de l'energie. (2021). Hydro [Dataset]. Retrieved from https://www.hydro.eaufrance.fr/rechercher/entites-hydrometriques

Ministerio para la Transición Ecológica y el Reto Demográfico. (2020). Geoportal [Dataset]. Retrieved from https://www.mitego.gob.es/es/cartografia-y-sig/ide/geoportal.html

Moriasi, D. N., Arnold, J. G., Liew, M. W. V., Bingner, R. L., Harmel, R. D., & Veith, T. L. (2007). Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. *Transactions of the ASABE*, **50**, 16.

Nguyen-Quang, T., Polcher, J., Ducharne, A., Arsouze, T., Zhou, X., Schneider, A., & Fita, L. (2018). ORCHIDEE-ROUTING: Revising the river routing scheme using a high-resolution hydrological database. *Geoscientific Model Development Discussions*, **11**(12), 4965–4985. https://doi.org/10.5194/gmd-11-4965-2018

Ning, T., Zhou, S., Chang, F., Shen, H., Li, Z., & Liu, W. (2019). Interaction of vegetation, climate and topography on evapotranspiration modelling at different time scales within the Budyko framework. *Agricultural and Forest Meteorology*, **275**, 59–68. https://doi.org/10.1016/j.agrformet.2019.05.001

Padrón, R. S., Gudmundsson, L., Greve, P., & Seneviratne, S. I. (2017). Large-scale controls of the surface water balance over land: Insights from a systematic review and meta-analysis. *Water Resources Research*, **53**(11), 9659–9678. https://doi.org/10.1002/2017WR021215

Polcher, J., Schrapffer, A., Dupont, E., Rinchiuso, L., Zhou, X., Boucher, O., et al. (2022). Hydrological modelling on atmospheric grids; using graphs of sub-grid elements to transport energy and water. *EGUSphere*, 1–34. https://doi.org/10.5194/egusphere-2022-690
Quintana-Seguí, P., Barella-Ortiz, A., Regueiro-Sanfiz, S., & Miguez-Macho, G. (2020). The utility of land-surface model simulations to provide drought information in a water management context using global and local forcing datasets. *Water Resources Management*, 34(7), 2135–2156.
https://doi.org/10.1007/s11269-018-2160-9

Quintana-Seguí, P., Turco, M., Herrera, S., & Miguez-Macho, G. (2017). Validation of a new SAFRAN-based gridded precipitation product for Spain and comparisons to Spain02 and ERA-Interim. *Hydrology and Earth System Sciences*, 21(4), 2187–2201. https://doi.org/10.5194/hess-21-2187-2017

Reaver, N. G. F., Kaplan, D. A., Klammler, H., & Jawitz, J. W. (2020). Technical note: Analytical inversion of the parametric Budyko equations. *Hydrology and Earth System Sciences Discussions*, 1–19. https://doi.org/10.5194/hess-2020-585

Reaver, N. G. F., Kaplan, D. A., Klammler, H., & Jawitz, J. W. (2022). Theoretical and empirical evidence against the Budyko catchment trajec-
tory conjecture. *Hydrology and Earth System Sciences*, 26(5), 1507–1525. https://doi.org/10.5194/hess-26-1507-2022

Ribes, A., Thao, S., Vautard, R., Dubuisson, B., Somot, S., Colin, J., et al. (2019). Observed increase in extreme daily rainfall in the French Mediterranean. *Climate Dynamics*, 52(1), 1095–1114. https://doi.org/10.1007/s00382-018-4179-2

Riedel, T., & Weber, T. K. D. (2020). Review: The influence of global change on Europe's water cycle and groundwater recharge. *Hydrogeology Journal*, 28(6), 1939–1959. https://doi.org/10.1007/s10040-020-02165-3

Roderick, M. L., & Farquhar, G. D. (2011). A simple framework for relating variations in runoff to variations in climatic conditions and catchment properties. *Water Resources Research*, 47(12). https://doi.org/10.1029/2010WR009826

Rottler, E., Francke, T., Bürger, G., & Bronstert, A. (2020). Long-term changes in central European river discharge for 1869–2016: Impact of changing snow covers, reservoir constructions and an intensified hydrological cycle. *Hydrology and Earth System Sciences*, 24(4), 1721–1740.
https://doi.org/10.5194/hess-24-1721-2020

Schneider, C., Laizé, C. L. R., Acreman, M. C., & Flörke, M. (2013). How will climate change modify river flow regimes in Europe? *Hydrology and Earth System Sciences*, 17(1), 325–339. https://doi.org/10.5194/hess-17-325-2013

Simons, G. W. H., Bastiaanssen, W. G. M., Cheema, M. J. M., Ahmad, B., & Immerzeel, W. W. (2020). A novel method to quantify consumed fractions and non-consumptive use of irrigation water: Application to the Indus Basin Irrigation System of Pakistan. *Agricultural Water Management*, 236, 106174. https://doi.org/10.1016/j.agwat.2020.106174

Stahl, K., Hisdal, H., Hannaford, J., Tallaksen, L. M., van Lanen, H. A. J., Sauquet, E., et al. (2010). Streamflow trends in Europe: Evidence from a dataset of near-natural catchments. *Hydrology and Earth System Sciences*, 14(12), 2367–2382. https://doi.org/10.5194/hess-14-2367-2010

Tafasca, S., Ducharne, A., & Valentin, C. (2020). Weak sensitivity of the terrestrial water budget to global soil texture maps in the ORCHIDEE land surface model. *Hydrology and Earth System Sciences*, 24(7), 3753–3774. https://doi.org/10.5194/hess-24-3753-2020

Tian, L., Jin, J., Wu, P., & Niu, G.-y. (2018). Quantifying the impact of climate change and human activities on streamflow in a semi-arid watershed with the Budyko equation incorporating dynamic vegetation information. *Water*, 10(12), 1781. https://doi.org/10.3390/w10121781

Tuel, A., Schaeфli, B., Zscheischler, J., & Martius, O. (2022). On the links between sub-seasonal clustering of extreme precipitation and high discharge in Switzerland and Europe. *Hydrology and Earth System Sciences*, 26(10), 2649–2669. https://doi.org/10.5194/hess-26-2649-2022

Vicente-Serrano, S. M., Lopez-Moreno, J.-I., Beguería, S., Lorenzo-Lacruz, J., Sanchez-Lorenzo, A., García-Ruiz, J. M., et al. (2014). Evidence of increasing drought severity caused by temperature rise in southern Europe. *Environmental Research Letters*, 9(4), 044001. https://doi.org/10.1088/1748-9326/9/4/044001

Vicente-Serrano, S. M., Peña-Gallardo, M., Hannaford, J., Murphy, C., Lorenzo-Lacruz, J., Dominguez-Castro, F., et al. (2019). Climate, irri-
gation, and land cover change explain streamflow trends in countries bordering the northeast Atlantic. *Geophysical Research Letters*, 46(19),
10821–10833. https://doi.org/10.1029/2019GL084084

Wang, F., Polcher, J., Peylin, P., & Bastrikov, V. (2018). Assimilation of river discharge in a land surface model to improve estimates of the conti-
nental water cycles. *Hydrology and Earth System Sciences*, 22(7), 3863–3882. https://doi.org/10.5194/hess-22-3863-2018

Wang, W., Zhang, Y., & Tang, Q. (2020). Impact assessment of climate change and human activities on streamflow signatures in the Yellow River Basin using the Budyko hypothesis and derived differential equation. *Journal of Hydrology*, 591, 125460. https://doi.org/10.1016/j.jhydrol.2020.125460

Weedon, G. P., Balsamo, G., Bellouin, N., Gomes, S., Best, M. J., & Viterbo, P. (2014). The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanalysis data. *Water Resources Research*, 50(9), 7505–7514.
https://doi.org/10.1002/2014WR015638

Xing, W., Wang, W., Shao, Q., & Yong, B. (2018). Identification of dominant interactions between climatic seasonality, catchment characteristics and agricultural activities on Budyko-type equation parameter estimation. *Journal of Hydrology*, 556, 585–599. https://doi.org/10.1016/j.jhydrol.2017.11.048

Xiong, M., Huang, C.-S., & Yang, T. (2020). Assessing the impacts of climate change and land use/cover change on runoff based on improved Budyko framework models considering arbitrary partition of the impacts. *Water*, 12(6), 1612. https://doi.org/10.3390/w12061612

Yang, D., Sun, F., Liu, Z., Cong, Z., Ni, G., & Lei, Z. (2007). Analyzing spatial and temporal variability of annual water-energy balance in nonhu-
mid regions of China using the Budyko hypothesis. *Water Resources Research*, 43(4). https://doi.org/10.1029/2006WR005224

Yang, H., Yang, D., Lei, Z., & Sun, F. (2008). New analytical derivation of the mean annual water-energy balance equation. *Water Resources Research*, 44(3). https://doi.org/10.1029/2007WR006135

Yang, Y., Zhang, S., McVicar, T. R., Beck, H. E., Zhang, Y., & Liu, B. (2018). Disconnection between trends of atmospheric drying and conti-
nental runoff. *Water Resources Research*, 54(7), 4700–4713. https://doi.org/10.1029/2018WR022593

Zhang, L., Hickel, K., Dawes, W. R., Chiew, F. H. S., Western, A. W., & Briggs, P. R. (2004). A rational function approach for estimating mean annual evapotranspiration. *Water Resources Research*, 40(2). https://doi.org/10.1029/2003WR002710

Zhang, L., Potter, N., Hickel, K., Zhang, Y., & Shao, Q. (2008). Water balance modeling over variable time scales based on the Budyko framework—Model development and testing. *Journal of Hydrology*, 360(1), 117–131. https://doi.org/10.1016/j.jhydrol.2008.07.021

Zhang, X., Dong, Q., Costa, V., & Wang, X. (2019). A hierarchical Bayesian model for decomposing the impacts of human activities and climate change on water resources in China. *Science of the Total Environment*, 665, 836–847. https://doi.org/10.1016/j.scitotenv.2019.02.189

Zhao, J., Huang, S., Huang, Q., Wang, H., & Leng, G. (2018). Detecting the dominant cause of streamflow decline in the Loess Plateau of China based on the latest Budyko equation. *Water*, 10(9), 1277. https://doi.org/10.3390/w10091277

Zheng, Y., Huang, Y., Zhou, S., Wang, K., & Wang, G. (2018). Effect partition of climate and catchment changes on runoff variation at the headwater region of the Yellow River based on the Budyko complementary relationship. *Science of the Total Environment*, 643, 1166–1177.
https://doi.org/10.1016/j.scitotenv.2018.06.195

Zveryaev, I. I. (2004). Seasonality in precipitation variability over Europe. *Journal of Geophysical Research: Atmospheres*, 109(D5), D05103.
https://doi.org/10.1029/2003JD003668