Research article 11 May 2020
Research article  11 May 2020
A lineintegralbased method to partition climate and catchment effects on runoff
 ^{1}NationalRegional Joint Engineering Research Center for Soil Pollution Control and Remediation in South China, Guangdong Key Laboratory of Integrated Agroenvironmental Pollution Control and Management, Guangdong Institute of Ecoenvironmental Science & Technology, Guangdong Academy of Sciences, Guangzhou 510650, China
 ^{2}Guangdong Engineering Center of Nonpoint Source Pollution Prevention Technology, Guangzhou 510650, China
 ^{1}NationalRegional Joint Engineering Research Center for Soil Pollution Control and Remediation in South China, Guangdong Key Laboratory of Integrated Agroenvironmental Pollution Control and Management, Guangdong Institute of Ecoenvironmental Science & Technology, Guangdong Academy of Sciences, Guangzhou 510650, China
 ^{2}Guangdong Engineering Center of Nonpoint Source Pollution Prevention Technology, Guangzhou 510650, China
Correspondence: Mingguo Zheng (mgzheng@soil.gd.cn)
Hide author detailsCorrespondence: Mingguo Zheng (mgzheng@soil.gd.cn)
It is a common task to partition the synergistic impacts of drivers in the environmental sciences. However, there is no mathematically precise solution to this partition task. Here I present a lineintegralbased method, which addresses the sensitivity to the drivers throughout the drivers' evolutionary paths so as to ensure a precise partition. The method reveals that the partition depends on both the change magnitude and pathway (timing of the change) but not on the magnitude alone unless used for a linear system. To illustrate this method, I applied the Budyko framework to partition the effects of climatic and catchment conditions on the temporal change in the runoff for 19 catchments from Australia and China. The proposed method reduces to the decomposition method when assuming a path in which climate change occurs first, followed by an abrupt change in catchment properties. The proposed method redefines the widely used sensitivity at a point as the pathaveraged sensitivity. The totaldifferential and the complementary methods simply concern the sensitivity at the initial and/or the terminal state, so they cannot give precise results. Although the pathaveraged sensitivities varied greatly among the catchments, they can be readily predicted within the Budyko framework. As a mathematically accurate solution, the proposed method provides a generic tool for conducting quantitative attribution analyses.
 Article
(1212 KB) 
Supplement
(320 KB)  BibTeX
 EndNote
The impacts of certain drivers on observed changes of interest often require quantification in environmental sciences. In the hydrology community, both climate and human activities have posed globalscale impact on hydrologic cycle and water resources (Barnett et al., 2008; Xu et al., 2014; Wang and Hejazi, 2011). Diagnosing their relative contributions to runoff is of considerable relevance to the researchers and managers. Unfortunately, performing a quantitative attribution analysis of runoff changes remains a challenge (Wang and Hejazi, 2011; Berghuijs and Woods, 2016; Zhang et al., 2016); this is to a considerable degree due to a lack of a mathematically precise method of decoupling synergistic and often confounding impacts of climate change and human activities.
Numerous studies have detected the longterm variability in runoff and attempted to partition the effects of climate change and human activities through various methods (Dey and Mishra, 2017); these include the pairedcatchments method and the hydrological modeling method. The pairedcatchment method can filter the effect of climatic variability and thus isolate the runoff change induced by vegetation changes (Brown et al., 2005). However, this method is capital intensive; moreover, it generally involves small catchments and experiences difficulties when extrapolating to large catchments (Zhang et al., 2011). Physically based hydrological models often have limitations such as a high data requirement, laborintensive calibration and validation processes, and inherent uncertainty and interdependence in parameter estimations (Binley et al., 1991; Wang et al., 2013; Liang et al., 2015). Conceptual models such as Budykotype equations (see Sect. 2.1) have consequently gained interest in recent years.
Within the Budyko framework, studies (Roderick and Farquhar, 2011; Zhang et al., 2016) have used the total differential of runoff as a proxy for the runoff change and the partial derivatives as the sensitivities (hereafter called the totaldifferential method). The total differential, however, is simply a firstorder approximation of the observed change (Fig. 1a). This approximation has caused an error in the calculation of climate impact on runoff, with the deviation ranging from 0 to $\mathrm{20}\times {\mathrm{10}}^{\mathrm{3}}$ m (or −118 to 174 %) in China (Yang et al., 2014). The elasticity method proposed by Schaake (1990) is also based on the totaldifferential expression (Sankarasubramanian et al., 2001; Zheng et al., 2009). The method uses the “elasticity” concept to assess the climate sensitivity of runoff. The elasticity coefficients, however, have been estimated in an empirical way and are not physically sound (Roderick and Farquhar, 2011; Liang et al., 2015).
The socalled decomposition method developed by Wang and Hejazi (2011) has also been widely used. The method assumes that climate changes cause a shift along a Budyko curve and then human interferences cause a vertical shift from one Budyko curve to another (Fig. 2). Under this assumption, the method extrapolates the Budyko models that are calibrated using observations of the reference period, in which human impacts remain minimal, to determine the humaninduced runoff changes that occur during the evaluation period.
Recently, Zhou et al. (2016) established a Budyko complementary relationship for runoff and further applied it to partitioning the climate and catchment effects. Superior to the totaldifferential method, the complementary method culminates by yielding a noresidual partition. Nevertheless, this method depends on a given weighted factor that is determined in an empirical but not a precise way. Furthermore, Zhou et al. (2016) argued that the partition is not unique in the Budyko framework because the path of the climate and catchment changes cannot be uniquely identified.
Obtaining a precise partition remains difficult, even when giving a precise mathematical model. This difficulty can be illustrated by using a precise hydrology model R=f (x, y), where R represents runoff, and x and y represent the climate factors and catchment characteristics, respectively. We assumed that R changes by ΔR when x changes by Δx and y changes by Δy; i.e., $\mathrm{\Delta}R=f(x+\mathrm{\Delta}x,y+\mathrm{\Delta}y)f(x,y)$. To determine the effect of x on ΔR – i.e., ΔR_{x} – a common practice is to assume that y remains constant when x changes by Δx. We thus obtain $\mathrm{\Delta}{R}_{x}=f(x+\mathrm{\Delta}x,y)f(x,y)$. Similarly, we can obtain $\mathrm{\Delta}{R}_{y}=f(x,y+\mathrm{\Delta}y)f(x,y)$. Although this derivation seems quite reasonable, it is problematic as $\mathrm{\Delta}{R}_{x}+\mathrm{\Delta}{R}_{y}\ne \mathrm{\Delta}R$. A further examination shows that a variable's effect on R seems to differ depending on the changing path (timing of the change). For example, $\mathrm{\Delta}{R}_{x}=f(x+\mathrm{\Delta}x,y)f(x,y)$ and $\mathrm{\Delta}{R}_{y}=f(x+\mathrm{\Delta}x,y+\mathrm{\Delta}y)f(x+\mathrm{\Delta}x,y)$ if x changes first and y subsequently changes (note that the partition is precise with $\mathrm{\Delta}{R}_{x}+\mathrm{\Delta}{R}_{y}=\mathrm{\Delta}R$ at this moment). If y changes first and x subsequently changes, the partition then becomes $\mathrm{\Delta}{R}_{x}=f(x+\mathrm{\Delta}x,y+\mathrm{\Delta}y)f(x,y+\mathrm{\Delta}y)$ and $\mathrm{\Delta}{R}_{y}=f(x,y+\mathrm{\Delta}y)f(x,y)$. In the case of x and y changing simultaneously, unfortunately, the current literature seems not to provide a mathematically precise solution.
The aim of this study is to propose a mathematically precise method to conduct a quantitative attribution to drivers. The method is based on the line integer (called the LI method hereafter) and takes into account the sensitivity throughout the evolutionary path of the drivers rather than at a point as the totaldifferential method does. To present and evaluate the proposed method, I decomposed the relative influences of climate and catchment conditions on runoff within the Budyko framework using data from 19 catchments from Australia and China.
2.1 Budyko framework and the MCY equation
Budyko (1974) argued that mean annual evapotranspiration (E) is largely determined by the water and energy balance of a catchment. Using precipitation (P) and potential evapotranspiration (E_{0}) as proxies for water and energy availabilities, respectively, the Budyko framework relates evapotranspiration losses to the aridity index defined as the ratio of E_{0} over P. The Budyko framework has gained wide acceptance in the hydrology community (Berghuijs and Woods, 2016; Sposito, 2017). In recent decades, several equations have been developed to describe the Budyko framework. Among them, the Mezentsev–Choudhury–Yang equation (Mezentsev, 1955; Choudhury, 1999; Yang et al., 2008) (called the MCY equation hereafter) has been widely accepted and was used in this study:
where $n\in (\mathrm{0},\mathrm{\infty})$ is an integration constant that is dimensionless and represents catchment properties. Equation (1) requires a relatively long timescale whereby the water storage of a catchment is negligible and the water balance equation reduces to be $R=PE$. Here I adopted a “tuned” n value that can obtain an exact accordance between the E calculated by Eq. (1) and that actually encountered (P−R).
The partial differentials of R with respect to P, E_{0}, and n are given as
2.2 Theory of the lineintegralbased method
We start by considering an example of a twovariable function $z=f(x,y)$ and assuming that x and y are independent. The function has continuous partial derivatives $\partial z/\partial x={f}_{x}(x,y)$ and $\partial z/\partial y={f}_{y}(x,y)$. Suppose that x and y vary along a smooth curve L (e.g., in Fig. 3) from the initial state (x_{0}, y_{0}) to the terminal state (x_{N}, y_{N}), and z covaries from z_{0} to z_{N}. Let $\mathrm{\Delta}z={z}_{N}{z}_{\mathrm{0}}$, $\mathrm{\Delta}x={x}_{N}{x}_{\mathrm{0}}$, and $\mathrm{\Delta}y={y}_{N}{y}_{\mathrm{0}}$. Our goal is to determine a mathematical solution that quantifies the effects of Δx and Δy on Δz, i.e., Δz_{x} and Δz_{y}. Δz_{x} and Δz_{y} should be subject to the constraint $\mathrm{\Delta}{z}_{x}+\mathrm{\Delta}{z}_{y}=\mathrm{\Delta}z$.
As shown in Fig. 3, points M_{1}(x_{1},y_{1}), …, ${M}_{N\mathrm{1}}({x}_{N\mathrm{1}},{y}_{N\mathrm{1}})$ partition L into N distinct segments. Let $\mathrm{\Delta}{x}_{i}={x}_{i+\mathrm{1}}{x}_{i}$, $\mathrm{\Delta}{y}_{i}={y}_{i+\mathrm{1}}{y}_{i}$, and $\mathrm{\Delta}{z}_{i}={z}_{i+\mathrm{1}}{z}_{i}$. For each segment, Δz_{i} can be approximated as dz_{i}: $\mathrm{\Delta}{z}_{i}\approx \mathrm{d}{z}_{i}={f}_{x}({x}_{i},{y}_{i})\mathrm{\Delta}{x}_{i}+{f}_{y}({x}_{i},{y}_{i})\mathrm{\Delta}{y}_{i}$. We then have $\mathrm{\Delta}z=\sum _{i=\mathrm{1}}^{N}\mathrm{\Delta}{z}_{i}\approx \sum _{i=\mathrm{1}}^{N}{f}_{x}({x}_{i},{y}_{i})\mathrm{\Delta}{x}_{i}+\sum _{i=\mathrm{1}}^{N}{f}_{y}({x}_{i},{y}_{i})\mathrm{\Delta}{y}_{i}$. We thus obtain the following respective approximation of Δz_{x} and Δz_{y}: $\mathrm{\Delta}{z}_{x}\approx \sum _{i=\mathrm{1}}^{N}{f}_{x}({x}_{i},{y}_{i})\mathrm{\Delta}{x}_{i}$ and $\mathrm{\Delta}{z}_{y}\approx \sum _{i=\mathrm{1}}^{N}{f}_{y}({x}_{i},{y}_{i})\mathrm{\Delta}{y}_{i}$. Next, define τ as the maximum length of the N segments. The smaller the value of τ, the closer the value of dz_{i} is to Δz_{i} and thus the more accurate the approximations are. The approximations become exact in the limit τ→0. Taking the limit τ→0 then converts the sum into integrals and gives a precise expression (this is an informal derivation; please see Appendix A for a formal one): $\mathrm{\Delta}z=\underset{\mathit{\tau}\to \mathrm{0}}{lim}\sum _{i=\mathrm{1}}^{N}{f}_{x}({x}_{i},{y}_{i})\mathrm{\Delta}{x}_{i}+\underset{\mathit{\tau}\to \mathrm{0}}{lim}\sum _{i=\mathrm{1}}^{N}{f}_{y}({x}_{i},{y}_{i})\mathrm{\Delta}{y}_{i}={\int}_{L}{f}_{x}(x,y)\mathrm{d}x+{\int}_{L}{f}_{y}(x,y)\mathrm{d}y$, where ${\int}_{L}{f}_{x}(x,y)\mathrm{d}x=\underset{\mathit{\tau}\to \mathrm{0}}{lim}\sum _{i=\mathrm{1}}^{N}{f}_{x}({x}_{i},{y}_{i})\mathrm{\Delta}{x}_{i}$ and ${\int}_{L}{f}_{y}(x,y)\mathrm{d}y=\underset{\mathit{\tau}\to \mathrm{0}}{lim}\sum _{i=\mathrm{1}}^{N}{f}_{y}({x}_{i},{y}_{i})\mathrm{\Delta}{y}_{i}$ denote the line integral of f_{x} and f_{y} along L (termed integral path) with respect to x and y, respectively. ${\int}_{L}{f}_{x}(x,y)\mathrm{d}x$ and ${\int}_{L}{f}_{y}(x,y)\mathrm{d}y$ exist provided that f_{x} and f_{y} are continuous along L. We thus obtain a precise evaluation of Δz_{x} and Δz_{y}:
Unlike the totaldifferential method, the sum of Δz_{x} and Δz_{y} persistently equals Δz (Appendix B). If f(x,y) is linear, then f_{x} and f_{y} are constant. Let f_{x}(x,y) and f_{y}(x,y) remain constant at C_{x} and C_{y}, respectively; then Δz_{x}=C_{x}Δx and Δz_{y}=C_{y}Δy. Δz_{x} and Δz_{y} are thus independent of L. If f(x,y) is nonlinear, however, both Δz_{x} and Δz_{y} vary with L, as is exemplified in Appendix C. Hence, the initial and the terminal states, together with the path connecting them, determine the resultant partition unless f(x,y) is linear.
The mathematical derivation above applies to a threevariable function as well. By doing the line integrals for the MCY equation, we obtain the desired results:
where ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} denote the effects on runoff change of P, E_{0}, and n, respectively. The sum of ΔR_{p} and $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$ represents the effect of climate change, and ΔR_{n} is often related to human activities although it probably includes the effects of other factors, such as climate seasonality (Roderick and Farquhar, 2011; Berghuijs and Woods, 2016). L denotes a threedimensional curve along which climate and catchment changes have occurred. I approximated L by a series of line segments. ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} were finally determined by summing up the integrals along each of the line segments (see Sect. 2.3).
2.3 Using the LI method to determine ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} within the Budyko framework

Determining ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} assuming a linear integral path:
A curve can always be approximated as a series of line segments. Hence, we can first handle the case of a linear integral path. Given two consecutive periods and assuming that the catchment state has evolved from (P_{1}, E_{01}, n_{1}) to (P_{2}, E_{02}, n_{2}) along a straight line L, let $\mathrm{\Delta}P={P}_{\mathrm{2}}{P}_{\mathrm{1}}$, $\mathrm{\Delta}{E}_{\mathrm{0}}={E}_{\mathrm{02}}{E}_{\mathrm{01}}$, and $\mathrm{\Delta}n={n}_{\mathrm{2}}{n}_{\mathrm{1}}$; then the line L is given by parametric equations: $P=\mathrm{\Delta}Pt+{P}_{\mathrm{1}}$, ${E}_{\mathrm{0}}=\mathrm{\Delta}{E}_{\mathrm{0}}t+{E}_{\mathrm{01}}$, $n=\mathrm{\Delta}nt+{n}_{\mathrm{1}}$, and $t\in [\mathrm{0},\mathrm{1}]$. Given these equations, Eq. (2) becomes a univariate function of t, i.e., $\partial R/\partial P={R}_{p}\left(t\right)$, $\partial R/\partial {E}_{\mathrm{0}}={R}_{{E}_{\mathrm{0}}}\left(t\right)$, and $\partial R/\partial n={R}_{n}\left(t\right)$. Then, ΔR_{P}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} can be evaluated as
$$\begin{array}{}\text{(5a)}& {\displaystyle}\begin{array}{rl}\mathrm{\Delta}{R}_{P}& ={\int}_{L}{\displaystyle \frac{\partial R}{\partial P}}\mathrm{d}P={\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{P}\left(t\right)\mathrm{d}(\mathrm{\Delta}Pt+{P}_{\mathrm{1}})\\ & =\mathrm{\Delta}P{\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{P}\left(t\right)\mathrm{d}t,\end{array}\text{(5b)}& {\displaystyle}\begin{array}{rl}\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}=& {\int}_{L}{\displaystyle \frac{\partial R}{\partial {E}_{\mathrm{0}}}}\mathrm{d}{E}_{\mathrm{0}}={\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{{E}_{\mathrm{0}}}\left(t\right)\mathrm{d}(\mathrm{\Delta}{E}_{\mathrm{0}}t+{E}_{\mathrm{01}})\\ & =\mathrm{\Delta}{E}_{\mathrm{0}}{\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{{E}_{\mathrm{0}}}\left(t\right)\mathrm{d}t,\end{array}\text{(5c)}& {\displaystyle}\begin{array}{rl}\mathrm{\Delta}{R}_{n}& ={\int}_{L}{\displaystyle \frac{\partial R}{\partial n}}\mathrm{d}n={\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{n}\left(t\right)\mathrm{d}(\mathrm{\Delta}nt+{n}_{\mathrm{1}})\\ & =\mathrm{\Delta}n{\int}_{\mathrm{0}}^{\mathrm{1}}{R}_{n}\left(t\right)\mathrm{d}t.\end{array}\end{array}$$Unfortunately, I could not determine the antiderivatives of R_{P}(t)dt, ${R}_{{E}_{\mathrm{0}}}\left(t\right)\mathrm{d}t$, and R_{n}(t)dt and had to make approximate calculations. As the discrete equivalent of integration is a summation, we can approximate the integration as a summation. I divided the $t\in [\mathrm{0},\mathrm{1}]$ interval into 1000 subintervals of the same width, i.e., setting dt identically equal to 0.001, and then calculated R_{P}(t)dt, ${R}_{{E}_{\mathrm{0}}}\left(t\right)\mathrm{d}t$, and R_{n}(t)dt for each subinterval. Let t_{i}=0.001i; $i\in [\mathrm{0},\mathrm{999}]$ and is integervalued. ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} are approximated as
$$\begin{array}{}\text{(6a)}& {\displaystyle}\mathrm{\Delta}{R}_{P}\approx \mathrm{0.001}\mathrm{\Delta}P\sum _{i=\mathrm{0}}^{\mathrm{999}}{R}_{P}\left({t}_{i}\right),\text{(6b)}& {\displaystyle}\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}\approx \mathrm{0.001}\mathrm{\Delta}{E}_{\mathrm{0}}\sum _{i=\mathrm{0}}^{\mathrm{999}}{R}_{{E}_{\mathrm{0}}}\left({t}_{i}\right),\text{(6c)}& {\displaystyle}\mathrm{\Delta}{R}_{n}\approx \mathrm{0.001}\mathrm{\Delta}n\sum _{i=\mathrm{0}}^{\mathrm{999}}{R}_{n}\left({t}_{i}\right).\end{array}$$ 
Dividing the evaluation period into a number of subperiods:
I first determined a change point and divided the whole observation period into the reference and evaluation periods. To determine the integral path, the evaluation period was further divided into a number of subperiods. The Budyko framework assumes a steadystate condition of a catchment and therefore requires no change in soil water storage. Over a time period of 5–10 years, it is reasonable to assume that changes in soil water storage will be sufficiently small (Zhang et al., 2001). Here, I divided the evaluation period into a number of 7year subperiods with the exception for the final subperiod, which varied from 7 to 13 years in length depending on the length of the evaluation period.

Determining ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} by approximating the integral path as a series of line segments:
For a short period, the integral path L can be considered as linear, which implies a temporally invariant change rate. For a long period in which the change rate varies over time, L can be fitted using a number of line segments. Given a reference period and an evaluation period comprising N subperiods, the catchment state was assumed to evolve from (P_{0}, E_{00}, n_{0}), …, (P_{i}, E_{0i}, n_{i}), …, to (P_{N}, E_{0N}, n_{N}), where the subscript “0” denotes the reference period, and “i” and “N” denote the ith and the final subperiods of the evaluation period, respectively. I used a series of line segments L_{1}, L_{2}, …, L_{N} to approximate the integral path L, where L_{i} connects points (P_{i−1}, ${E}_{\mathrm{0},i\mathrm{1}}$, n_{i−1}) with (P_{i}, E_{0i}, n_{i}). Then ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} were evaluated as the sum of the integrals along each of the line segments, which were calculated using Eq. (6).
^{a} R, P, and E_{0} represent the mean annual runoff, precipitation, and potential evaporation, all in 10^{−3}m yr^{−1}. n (dimensionless) is the parameter representing catchment properties in the MCY equation. AI is the dimensionless aridity index ($\mathit{\text{AI}}={E}_{\mathrm{0}}/P$). Data of catchments 1–12 were derived from Zhang et al. (2010). Data of catchments 13–16 were from Sun et al. (2014). Data of catchments 17, 18, and 19 were from Zheng et al. (2009), Jiang et al. (2015), and Gao et al. (2016), respectively. I used the change points given in the literatures to divide the observation period into the reference and elevation periods. The LI method further divides the evaluation period into a number of subperiods. The column “The final subperiod” denotes the final subperiod, which was used as the evaluation period for the totaldifferential method, the decomposition method, and the complementary method. ^{b} Catchments 1–12 are in Australia, and the others are in China. (1) Adjungbilly Creek; (2) Batalling Creek; (3) Bombala River; (4) Crawford River; (5) Darlot Creek; (6) Eumeralla River; (7) Goobarragandra Creek; (8) Jingellic Creek; (9) Mosquito Creek; (10) Traralgon Creek; (11) upper Denmark River; (12) Yate Flat Creek; (13) Yangxian station, Han River; (14) Ankang station, Han River; (15) Baihe station, Han River; (16) Danjiangkou station, Han River; (17) headwaters of the Yellow River basin; (18) Wei River; (19) Yan River.
2.4 Totaldifferential, decomposition, and complementary methods
To evaluate the LI method, I compared it with the existing methods, including the decomposition method, the totaldifferential method, and the complementary method. The totaldifferential method approximated ΔR as dR:
where ${\mathit{\lambda}}_{P}=\partial R/\partial P$, ${\mathit{\lambda}}_{{E}_{\mathrm{0}}}=\partial R/\partial {E}_{\mathrm{0}}$, and ${\mathit{\lambda}}_{n}=\partial R/\partial n$, representing the sensitivity coefficient of R with respect to P, E_{0}, and n, respectively. Within the totaldifferential method,ΔR_{P}=λ_{P}ΔP, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}={\mathit{\lambda}}_{{E}_{\mathrm{0}}}\mathrm{\Delta}{E}_{\mathrm{0}}$, and ΔR_{n}=λ_{n}Δn. I used the forward approximation, i.e., substituting the observed mean annual values of the reference period into Eq. (2), to estimate λ_{P}, ${\mathit{\lambda}}_{{E}_{\mathrm{0}}}$, and λ_{n}, as is standard in most studies (Roderick and Farquhar, 2011; Yang and Yang, 2011; Sun et al., 2014).
The decomposition method (Wang and Hejazi, 2011) calculated ΔR_{n} as follows:
where R_{2}, P_{2}, and E_{2} represent the mean annual runoff, precipitation, and evapotranspiration of the evaluation period, respectively; R_{2}^{′} and E_{2}^{′} represent the mean annual runoff and evapotranspiration, respectively, given the climate conditions of the evaluation period and the catchment conditions of the reference period (Fig. 2). Both E_{2} and E_{2}^{′} were calculated by Eq. (1) but using n values of the evaluation period and the reference period, respectively.
The complementary method (Zhou et al., 2016) uses a linear combination of the complementary relationship for runoff to determine ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n}:
where the subscripts “1” and “2” denotes the reference and the evaluation periods, respectively. a is a weighting factor and varies from 0 to 1. As suggested by Zhou et al. (2016), I set a=0.5. Equation (9) thus gave an estimation of ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} as follows:
2.5 Data
I collected runoff and climate data from 19 selected catchments evaluated in previous studies (Table 1). The changepoint years given in these studies were directly used to determine the reference and evaluation periods for the LI method. As mentioned above, the LI method further divides the evaluation period into a number of subperiods. For the sake of comparison, the final subperiod of the evaluation period was used as the evaluation period for the decomposition, the totaldifferential, and the complementary methods (it can be equally considered that all of the four methods used the final subperiod as the evaluation period, but the LI method cares about the intermediate period between the reference and the evaluation periods, and the other methods do not). Eight of the 19 catchments had a reference period comprising only one subperiod (Table 1), and the others had two to seven subperiods.
The 19 selected catchments have diverse climates and landscapes, with 12 from Australia and seven from China (Table 1). The catchments span from tropical to subtropical and temperate areas and from humid to semihumid and semiarid regions, with the mean annual rainfall varying from 506 to $\mathrm{1014}\times {\mathrm{10}}^{\mathrm{3}}$ m and potential evaporation from 768 to $\mathrm{1169}\times {\mathrm{10}}^{\mathrm{3}}$ m. The dryness index ranges between 0.86 and 1.91. The catchment areas vary by 5 orders of magnitude from 1.95 to 121 972 with a median 606×10^{6} m^{2}. The key data include annual runoff, precipitation, and potential evaporation. The record length varied between 19 and 76 with a median of 39 years. All the catchments experienced changes in climate and catchment properties over the observation periods. The precipitation changes from the reference to the evaluation period ranged between −153 and $\mathrm{79}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1}, and between −35 and $\mathrm{41}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1} for potential evaporation (Table 2). The coeval change in the parameter n of the MCY equation ranged between −0.2 and 1.4. The mean annual streamflow reduced for all catchments, ranging from 0.4 to 169 with a median $\mathrm{38}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1}. The change in catchment properties mainly refers to the vegetation cover or land use change. More details on the data and the catchments can be found in Zhang et al. (2010, 2011), Sun et al. (2014), Zheng et al. (2009), Jiang et al. (2015), and Gao et al. (2016).
Table 3 lists the resultant values of ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} from the LI method and the three other methods. Please see the Supplement for detailed calculation steps.
Figure 4a compares the resultant ΔR_{n} of the LI method and the decomposition method. Although they are quite similar, the discrepancies can be $\mathit{>}\mathrm{20}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1}. The decomposition method assumes that climate change occurs first and then human interferences cause a sudden change in catchment properties (Fig. 2). Such a fictitious path is identical to the path ABC in Fig. 3, provided that x represents climate factors and y catchment properties. When ABC was adopted as the integral path, the LI method yielded the same results as the decomposition method did (Fig. 4b). Hence, the decomposition method can be considered as a case of the LI method that uses a special integral path.
The totaldifferential method is predicated on an approximate equation, i.e., Eq. (7). The LI method reveals that the precise form of the equation is $\mathrm{\Delta}R=\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}\mathrm{\Delta}P+\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}\mathrm{\Delta}{E}_{\mathrm{0}}+\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}\mathrm{\Delta}n$ (Appendix D), where $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ denote the pathaveraged sensitivity of R to P, E_{0}, and n, respectively. All points along the path have the same weight in determining $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$. To determine them, the totaldifferential and the complementary methods utilize only the initial and/or the terminal states. Neglecting the intermediate states results in an imprecise partition, as was illustrated in Figs. 5 and 6 as well as in Fig. 1 using a univariate function, and even a reverse trend estimation (see $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$ for catchment no. 1 in Table 3).
As with the LI method, the complementary method produced ΔR_{P}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} that exactly summed up to ΔR. Although its resultant ΔR_{p}, $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}$, and ΔR_{n} values were all in accordance with the LI method (Fig. 6), the LI method often yielded values beyond the bounds given by the complementary method (Fig. 7); this is because the maximum or minimum sensitivities do not necessarily occur at the initial or terminal states.
$\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}},\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ imply the average runoff change induced by a unit change in P, E_{0}, and n, respectively (Appendix D). $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ all varied by up to several times or even 10fold between the studied catchments (Table S4). Figure 8 shows that Eq. (2) reproduced $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ very well when taking the longterm means of P, E_{0}, and n as inputs, a reflection of the fact that the dependent variable approached its average when the independent variables were set to be their averages. This finding is of relevance to the spatial prediction of $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$.
The LI method highlights the role of the evolutionary path in determining the resultant partition. Yet it seems that no studies have accounted for the path issue while evaluating the relative influences of drivers. The limit of the LI method is high data requirements for obtaining the evolutionary path. When the path data are unavailable, the complementary method can be considered as an alternative. The complementary method is free of residuals; moreover, it employs data from both the reference and the evaluation periods, thereby generally yielding sensitivities closer to the pathaveraged results than the totaldifferential method.
While using the Budyko models, a reasonable timescale is relevant to meet the assumption that changes in catchment water storage are small relative to the magnitude of fluxes of P, R, and E (Donohue et al., 2007; Roderick and Farquhar, 2011). A 7year timescale was used in the present study, as most studies have suggested that a time period of 5–10 years (Zhang et al., 2001, 2016; Wu et al., 2017a, b; Li et al., 2017) or even 1 year (Roderick and Farquhar, 2011; Sivapalan et al., 2011; Carmona et al., 2014; Ning et al., 2017) is reasonable. Nevertheless, some studies have argued that the time period should be longer than 10 years (Li et al., 2016; Dey and Mishra, 2017). Using the Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry, Zhao et al. (2011) detected the water storage variations for the three largest river basins of China, namely, the Yellow, Yangtze, and Zhujiang. The Yellow River mostly drains an arid and semiarid region (P, $\mathrm{450}\times {\mathrm{10}}^{\mathrm{3}}$ m; R, $\mathrm{70}\times {\mathrm{10}}^{\mathrm{3}}$ m; E, $\mathrm{380}\times {\mathrm{10}}^{\mathrm{3}}$ m), and the Yangtze (P, $\mathrm{110}\times {\mathrm{10}}^{\mathrm{3}}$ m; R, $\mathrm{550}\times {\mathrm{10}}^{\mathrm{3}}$ m; E, $\mathrm{550}\times {\mathrm{10}}^{\mathrm{3}}$ m) and the Zhujiang River basins (P, $\mathrm{1400}\times {\mathrm{10}}^{\mathrm{3}}$ m; R, $\mathrm{780}\times {\mathrm{10}}^{\mathrm{3}}$ m; E, $\mathrm{620}\times {\mathrm{10}}^{\mathrm{3}}$ m) are humid. The amplitude of the water storage variations between years was 7, 37.2, and $\mathrm{65}\times {\mathrm{10}}^{\mathrm{3}}$ m for the three rivers, respectively, at 1 order of magnitude smaller than the fluxes of P, R and E. Although the observations cannot be directly extrapolated to other regions, the possibility seems remote that the use of a 7year aggregated time period strongly violates the assumption of the steadystate condition.
The mutual independence of the drivers is crucial for a valid partition. In the present study, although annual P and E_{0} exhibited significant correlation for most catchments (p<0.05), the aggregated P, E_{0}, and n over a 7year period showed minimal correlation (mostly p>0.1). The interdependence between the drivers can considerably confound the resultant partitions of the LI method and other existing methods.
The LI method revises the concept of sensitivity at a point as the pathaveraged sensitivity. Mathematically, the LI method is unrelated to a functional form and hence applies to communities other than just hydrology. For example, identifying the carbon emission budgets (an allowable amount of anthropogenic CO_{2} emission consistent with a limiting warming target) is crucial for global efforts to mitigate climate change. The LI method suggested that the emission budgets depends on both the emission magnitude and pathway (timing of emissions), which is in line with a recent study by Gasser et al. (2018). An optimal pathway would facilitate an elevated carbon budget unless the carbon–climate system behaves in a linear fashion.
This study presented the LI method using timeseries data, but it applies equally to the case of spatial series of data. Given a model that relates fluvial or eolian sediment load to the influencing factors (e.g., rainfall and topography), for example, the LI method can be used to separate their contributions to the sedimentload change along a river or in the alongwind direction.
Based on the line integral, I created a mathematically precise method to partition the synergistic effects of several factors that cumulatively drive a system to change from one state to the other. The method is relevant for quantitative assessments of the relative roles of the factors behind the change in the system state. I applied the LI method to partition the effects of climatic and catchment conditions on runoff within the Budyko framework. The method reveals that, in addition to the change magnitude, the change pathways of climatic and catchment conditions also play a role. Instead of using the runoff sensitivity at a point, the LI method uses the pathaveraged sensitivity, thereby ensuring a mathematically precise partition. As a mathematically accurate scheme, the LI method has the potential to be a generic attribution approach in the environmental sciences.
Proof. In order to prove the equation above, we define the curve L in Fig. 3 as being given by a parametric equation: x=x(t), y=y(t), $t\in [{t}_{\mathrm{0}},{t}_{N}]$, and therefore $\mathrm{\Delta}z={z}_{N}{z}_{\mathrm{0}}=f\left[x\right({t}_{N}),y({t}_{N}\left)\right]f\left[x\right({t}_{\mathrm{0}}),y({t}_{\mathrm{0}}\left)\right]$. Substituting the parametric
equations, the righthand side of the equation can be rewritten as
${\int}_{L}{f}_{x}(x,y)\mathrm{d}x+{\int}_{L}{f}_{y}(x,y)\mathrm{d}y$
Let $g\left(t\right)=f\left[x\right(t),y(t\left)\right]$, and after using the chain rule to differentiate g with respect to t we obtain
Thus, g^{′}(t) is just the integrand in Eq. (A1), and the righthand side of the equation can be further rewritten as follows.
The righthand side of the equation $={\int}_{{t}_{\mathrm{0}}}^{{t}_{N}}{g}^{\prime}\left(t\right)\mathrm{d}t={\left[g\left(t\right)\right]}_{{t}_{\mathrm{0}}}^{{t}_{N}}=g\left({t}_{N}\right)g\left({t}_{\mathrm{0}}\right)=f\left[x\right({t}_{N}),y({t}_{N}\left)\right]f\left[x\right({t}_{\mathrm{0}}),y({t}_{\mathrm{0}}\left)\right]=$ the lefthand side of the equation. Q.E.D.
Preliminary theorem: Given an open, simply connected region G (i.e., no holes in G) and two functions P(x,y) and Q(x,y) that have continuous firstorder derivatives, if and only if $\partial P/\partial y=\partial Q/\partial x$ throughout G, then ${\int}_{L}P(x,y)\mathrm{d}x+{\int}_{L}Q(x,y)\mathrm{d}y$ is path independent; i.e., it depends solely on the starting and ending points of L.
Proof. We have $\partial {f}_{x}/\partial y={\partial}^{\mathrm{2}}z/\partial x\partial y$ and $\partial {f}_{y}/\partial x={\partial}^{\mathrm{2}}z/\partial y\partial x$. As ${\partial}^{\mathrm{2}}z/\partial x\partial y={\partial}^{\mathrm{2}}z/\partial y\partial x$, we can state that $\partial {f}_{x}/\partial y=\partial {f}_{y}/\partial x$, meeting the above condition and proving that ${\int}_{L}{f}_{x}(x,y)\mathrm{d}x+{\int}_{L}{f}_{y}(x,y)\mathrm{d}y$ is path independent. The statement was further exemplified using a fictitious example in Appendix C. Q.E.D.
Runoff (R, 10^{−3} m yr^{−1}) at a site is assumed to increase from 120 to $\mathrm{195}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1} with $\mathrm{\Delta}R=\mathrm{75}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1}; meanwhile, precipitation (P, 10^{−3} m yr^{−1}) varies from 600 to $\mathrm{650}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1} ($\mathrm{\Delta}P=\mathrm{75}\times {\mathrm{10}}^{\mathrm{3}}$ m yr^{−1}) and the runoff coefficient (C_{R}, dimensionless) varies from 0.2 to 0.3 (ΔC_{R}=0.1). The goal is to partition ΔR into the effects of the precipitation (ΔR_{P}) and runoff coefficient ($\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$), provided that P and C_{R} are independent. We have a function R=PC_{R} and its partial derivatives $\partial R/\partial P={C}_{\mathrm{R}}$ and $\partial R/\partial {C}_{\mathrm{R}}=P$. Given a path L along which P and C_{R} change and using Eq. (3), the LI method evaluates ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ as
The result differs depending on L, but the sum of ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ uniformly equals ΔR. This dynamic is demonstrated using Fig. 3, in which we considered that the x axis represents C_{R}, and the y axis P. Point A denotes the initial state (C_{R}=0.2, P=600), and point C the terminal state (C_{R}=0.3, P=650). I calculated ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ along three fictitious paths as follows:

L= AC. Line segment AC has equation $P=\mathrm{500}{C}_{\mathrm{R}}+\mathrm{500}$, $\mathrm{0.2}\le {C}_{\mathrm{R}}\le \mathrm{0.3}$. Let us take C_{R} as the parameter and write the equation in the parametric form as $P=\mathrm{500}{C}_{\mathrm{R}}+\mathrm{500}$, C_{R}=C_{R}, $\mathrm{0.2}\le {C}_{\mathrm{R}}\le \mathrm{0.3}$. By substituting the equation into Eq. (C1), we have
$$\begin{array}{rl}\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}& ={\int}_{\mathrm{AC}}P\mathrm{d}{C}_{\mathrm{R}}={\int}_{\mathrm{0.2}}^{\mathrm{0.3}}(\mathrm{500}{C}_{\mathrm{R}}+\mathrm{500})\mathrm{d}{C}_{\mathrm{R}}=\mathrm{62.5},\\ \mathrm{\Delta}{R}_{P}& ={\int}_{\mathrm{AC}}{C}_{\mathrm{R}}\mathrm{d}P={\int}_{\mathrm{AC}}{C}_{\mathrm{R}}\mathrm{d}(\mathrm{500}{C}_{\mathrm{R}}+\mathrm{500})\\ & =\mathrm{500}{\int}_{\mathrm{0.2}}^{\mathrm{0.3}}{C}_{\mathrm{R}}\mathrm{d}{C}_{\mathrm{R}}=\mathrm{12.5}.\end{array}$$ 
L= AB + BC. To evaluate on the broken line, we can evaluate separately on AB and BC and then sum them up. The equation for AB is P=600, $\mathrm{0.2}\le {C}_{\mathrm{R}}\le \mathrm{0.3}$, while for BC it is C_{R}=0.3, $\mathrm{600}\le P\le \mathrm{650}$. Note that a constant C_{R} or P implies that dC_{R}=0 or dP=0. Eq. (C1) then becomes
$$\begin{array}{rl}\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}& ={\int}_{\mathrm{AB}+\mathrm{BC}}P\mathrm{d}{C}_{\mathrm{R}}={\int}_{\mathrm{AB}}P\mathrm{d}{C}_{\mathrm{R}}+{\int}_{\mathrm{BC}}P\mathrm{d}{C}_{\mathrm{R}}\\ & ={\int}_{\mathrm{0.2}}^{\mathrm{0.3}}\mathrm{600}\mathrm{d}{C}_{\mathrm{R}}+\mathrm{0}=\mathrm{60},\\ \mathrm{\Delta}{R}_{P}& ={\int}_{\mathrm{AB}+\mathrm{BC}}{C}_{\mathrm{R}}\mathrm{d}P={\int}_{\mathrm{AB}}{C}_{\mathrm{R}}\mathrm{d}P+{\int}_{\mathrm{BC}}{C}_{\mathrm{R}}\mathrm{d}P\\ & =\mathrm{0}+{\int}_{\mathrm{600}}^{\mathrm{650}}\mathrm{0.3}\mathrm{d}P=\mathrm{15}.\end{array}$$ 
$L=\mathrm{AD}+\mathrm{DC}$. The equation for AD is C_{R}=0.2, $\mathrm{600}\le P\le \mathrm{650}$, and for DC it is P=650, $\mathrm{0.2}\le {C}_{\mathrm{R}}\le \mathrm{0.3}$. ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ are evaluated as
$$\begin{array}{rl}\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}& ={\int}_{\mathrm{AD}+\mathrm{DC}}P\mathrm{d}{C}_{\mathrm{R}}={\int}_{\mathrm{AD}}P\mathrm{d}{C}_{\mathrm{R}}+{\int}_{\mathrm{DC}}P\mathrm{d}{C}_{\mathrm{R}}\\ & =\mathrm{0}+{\int}_{\mathrm{0.2}}^{\mathrm{0.3}}\mathrm{650}\mathrm{d}{C}_{\mathrm{R}}=\mathrm{65},\\ \mathrm{\Delta}{R}_{P}& ={\int}_{\mathrm{AD}+\mathrm{DC}}{C}_{\mathrm{R}}\mathrm{d}P={\int}_{\mathrm{AD}}{C}_{\mathrm{R}}\mathrm{d}P+{\int}_{\mathrm{DC}}{C}_{\mathrm{R}}\mathrm{d}P\\ & ={\int}_{\mathrm{600}}^{\mathrm{650}}\mathrm{0.2}\mathrm{d}P+\mathrm{0}=\mathrm{10}.\end{array}$$As expected, the sum of ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ persistently equals ΔR although ΔR_{P} and $\mathrm{\Delta}{R}_{{C}_{\mathrm{R}}}$ varies with L.
If the interval [x_{0}, x_{N}] in Fig. 3 is partitioned into N distinct bins of the same width $\mathrm{\Delta}{x}_{i}=\mathrm{\Delta}x/N$, Eq (3a) can then be rewritten as
where $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{x}}=\underset{\mathit{\tau}\to \mathrm{0}}{lim}\frac{\sum _{i=\mathrm{1}}^{N}{f}_{x}({x}_{i},{y}_{i})}{N}$, denoting the average of f_{x}(x, y) along the curve L. Likewise, we have $\mathrm{\Delta}{Z}_{y}=\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{y}}\mathrm{\Delta}y$, where $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{y}}$ denotes the average of f_{y}(x,y) along the curve L. As a result
The result can readily be extended to a function of three variables. Applying the mathematic derivation determined above to the MCY equation results in a precise form of Eq. (7):
where $\mathrm{\Delta}{R}_{P}=\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}\mathrm{\Delta}P$; $\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}=\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}\mathrm{\Delta}{E}_{\mathrm{0}}$; $\mathrm{\Delta}{R}_{n}=\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}\mathrm{\Delta}n$; and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ denote the arithmetic mean of $\partial R/\partial P$, $\partial R/\partial {E}_{\mathrm{0}}$, and $\partial R/\partial n$ along a path of climate and catchment changes, respectively. Because $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}=\mathrm{\Delta}{R}_{P}/\mathrm{\Delta}P$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}=\mathrm{\Delta}{R}_{{E}_{\mathrm{0}}}/\mathrm{\Delta}{E}_{\mathrm{0}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}=\mathrm{\Delta}{R}_{n}/\mathrm{\Delta}n$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{P}}$, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{{E}_{\mathrm{0}}}}$, and $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{n}}$ also imply the runoff change due to a unit change in P, E_{0}, and n, respectively.
Given a onedimensional function z=f(x) and its derivative f^{′}(x), we assumed that f^{′}(x) averages $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{x}}$ over the range $(x,x+\mathrm{\Delta}x)$, i.e., $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{x}}=\underset{\mathit{\tau}\to \mathrm{0}}{lim}\frac{\sum _{i=\mathrm{1}}^{N}{f}^{\prime}\left({x}_{i}\right)}{N}$. According to the mean value theorem for integrals, $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{x}}={\int}_{x}^{x+\mathrm{\Delta}x}{f}^{\prime}\left(x\right)\mathrm{d}x/\mathrm{\Delta}x$. In terms of the Newton–Leibniz formula, ${\int}_{x}^{x+\mathrm{\Delta}x}{f}^{\prime}\left(x\right)\mathrm{d}x=f(x+\mathrm{\Delta}x)f\left(x\right)=\mathrm{\Delta}z$. Thus, we obtain $\stackrel{\mathrm{\u203e}}{{\mathit{\lambda}}_{x}}=\mathrm{\Delta}z/\mathrm{\Delta}x$.
The data used in this study are freely available by contacting the author.
The supplement related to this article is available online at: https://doi.org/10.5194/hess2423652020supplement.
MZ designed the study, analyzed the data, and wrote the manuscript.
The author declares that there is no conflict of interest.
I thank Yuanqiong Zheng for his assistance with the mathematical proof in Appendix A.
This research has been supported by the National Natural Science Foundation of China (grant no. 41671278) and the Guangdong Academy of Sciences' Project of Science and Technology Development (grant nos. 2019GDASYL0103043, 2019GDASYL0502004, 2019GDASYL0401003, 2019GDASYL0301002, 2019GDASYL0503003, and 2019GDASYL0102002).
This paper was edited by Erwin Zehe and reviewed by two anonymous referees.
Barnett, T. P., Pierce, D. W., Hidalgo, H. G., Bonfils, C., Santer, B. D., Das, T., .Bala G., Woods, A. W., Nozawa, T., Mirin, A. A., Cayan D. R., and Dettinger, M. D.: Humaninduced changes in the hydrology of the western United States, Science, 319, 1080–1083, https://doi.org/10.1126/science.1152538, 2008.
Berghuijs, W. R., and Woods, R. A.: Correspondence: Spacetime asymmetry undermines water yield assessment, Nat. Commun., 7, 11603, https://doi.org/10.1038/ncomms11603, 2016.
Binley, A. M., Beven, K. J., Calver, A., and Watts, L. G.: Changing responses in hydrology: assessing the uncertainty in physically based model predictions, Water Resour. Res., 27, 1253–1261, https://doi.org/10.1029/91WR00130, 1991.
Brown, A. E., Zhang, L., McMahon, T. A., Western, A. W., and Vertessy, R. A.: A review of paired catchment studies for determining changes in water yield resulting from alterations in vegetation, J. Hydrol., 310, 26–61, https://doi.org/10.1016/j.jhydrol.2004.12.010, 2005.
Budyko, M. I.: Climate and Life, Academic, NY, 508 pp., 1974.
Carmona, A. M., Sivapalan, M., Yaeger, M. A., and Poveda, G.: Regional patterns of interannual variability of catchment water balances across the continental US: A Budyko framework, Water Resour. Res., 50, 9177–9193, https://doi.org/10.1002/2014wr016013, 2014.
Choudhury, B. J.: Evaluation of an empirical equation for annual evaporation using field observations and results from a biophysical model, J. Hydrol., 216, 99–110, https://doi.org/10.1016/S00221694(98)002935, 1999.
Dey, P. and Mishra, A.: Separating the impacts of climate change and human activities on streamflow: A review of methodologies and critical assumptions, J. Hydrol., 548, 278–290, https://doi.org/10.1016/j.jhydrol.2017.03.014, 2017.
Donohue, R. J., Roderick, M. L., and McVicar, T. R.: On the importance of including vegetation dynamics in Budyko's hydrological model, Hydrol. Earth Syst. Sci., 11, 983–995, https://doi.org/10.5194/hess119832007, 2007.
Gao, G., Ma, Y., and Fu, B.: Multitemporal scale changes of streamflow and sediment load in a loess hilly watershed of China, Hydrol. Process., 30, 365–382, https://doi.org/10.1002/hyp.10585, 2016.
Gasser, T., Kechiar, M., Ciais, P., Burke, E. J., Kleinen, T., Zhu, D., Huang, Y., Ekici, A., and Obersteiner, M.: Pathdependent reductions in CO_{2} emission budgets caused by permafrost carbon release, Nat. Geosci., 11, 830–835, https://doi.org/10.1038/s4156101802270, 2018.
Jiang, C., Xiong, L., Wang, D., Liu, P., Guo, S., and Xu C.: Separating the impacts of climate change and human activities on runoff using the Budykotype equations with timevarying parameters, J. Hydrol., 522, 326–338, https://doi.org/10.1016/j.jhydrol.2014.12.060, 2015.
Li, Z., Ning, T., Li, J., and Yang, D.: Spatiotemporal variation in the attribution of streamflow changes in a catchment on China's Loess Plateau, Catena, 158, 1–8, https://doi.org/10.1016/j.catena.2017.06.008, 2017.
Liang, W., Bai, D., Wang, F., Fu, B., Yan, J., Wang, S., Yang, Y., Long, D., and Feng, M.: Quantifying the impacts of climate change and ecological restoration on streamflow changes based on a Budyko hydrological model in China's Loess Plateau, Water Resour. Res., 51, 6500–6519, https://doi.org/10.1002/2014WR016589, 2015.
Mezentsev, V. S.: More on the calculation of average total evaporation, Meteorol. Gidrol., 5, 24–26, 1955.
Ning, T., Li, Z., and Liu, W.: Vegetation dynamics and climate seasonality jointly control the interannual catchment water balance in the Loess Plateau under the Budyko framework, Hydrol. Earth Sys. Sci., 21, 1515–1526, https://doi.org/10.5194/hess2016484, 2017.
Roderick, M. L. and Farquhar, G. D.: A simple framework for relating variations in runoff to variations inclimatic conditions and catchment properties, Water Resour. Res., 47, W00G07, https://doi.org/10.1029/2010WR009826, 2011.
Sankarasubramanian, A., Vogel, R. M., and Limbrunner, J. F.: Climate elasticity of streamflow in the United States, Water Resour. Res., 37, 1771–1781, https://doi.org/10.1029/2000WR900330, 2001.
Schaake, J. C.: From climate to flow, Climate Change and U.S. Water Resources, edited by: Waggoner, P. E., chap. 8, John Wiley, NY, 177–206, 1990.
Sivapalan, M., Yaeger, M. A., Harman, C. J., Xu, X. Y., and Troch, P. A.: Functional model of water balance variability at the catchment scale: 1. Evidence of hydrologic similarity and spacetime symmetry, Water Resour. Res., 47, W02522, https://doi.org/10.1029/2010wr009568, 2011.
Sposito, G.: Understanding the Budyko equation, Water, 9, 236, https://doi.org/10.3390/w9040236, 2017.
Sun, Y., Tian, F., Yang, L., and Hu, H.: Exploring the spatial variability of contributions from climate variation and change in catchment properties to streamflow decrease in a mesoscale basin by three different methods, J. Hydrol., 508, 170–180, https://doi.org/10.1016/j.jhydrol.2013.11.004, 2014.
Wang, D. and Hejazi, M.: Quantifying the relative contribution of the climate and direct human impacts on mean annual streamflow in the contiguous United States, Water Resour. Res., 47, W00J12, https://doi.org/10.1029/2010WR010283, 2011.
Wang, W., Shao, Q., Yang, T., Peng, S., Xing, W., Sun, F., and Luo, Y.: Quantitative assessment of the impact of climate variability and human activities on runoff changes: A case study in four catchments of the Haihe River basin, China, Hydrol. Process., 27, 1158–1174, https://doi.org/10.1002/hyp.9299, 2013.
Wu, J., Miao, C., Wang, Y., Duan, Q., and Zhang, X.: Contribution analysis of the longterm changes in seasonal runoff on the Loess Plateau, China, using eight Budykobased methods, J. Hydrol., 545, 263–275, https://doi.org/10.1016/j.jhydrol.2016.12.050, 2017a.
Wu, J., Miao, C., Zhang, X., Yang, T., and Duan, Q.: Detecting the quantitative hydrological response to changes in climate and human activities, Sci. Total Environ., 586, 328–337, https://doi.org/10.1016/j.scitotenv.2017.02.010, 2017b.
Xu, X., Yang, D., Yang, H., and Lei, H.: Attribution analysis based on the Budyko hypothesis for detecting the dominant cause of runoff decline in Haihe basin, J. Hydrol., 510, 530–540, https://doi.org/10.1016/j.jhydrol.2013.12.052, 2014.
Yang, H. and Yang, D.: Derivation of climate elasticity of runoff to assess the effects of climate change on annual runoff, Water Resour. Res., 47, W07526, https://doi.org/10.1029/2010WR009287, 2011.
Yang, H., Yang, D., Lei, Z., and Sun, F.: New analytical derivation of the mean annual waterenergy balance equation, Water Resour. Res., 44, W03410, https://doi.org/10.1029/2007WR006135, 2008.
Yang, H., Yang, D., and Hu, Q.: An error analysis of the Budyko hypothesis for assessing the contribution of climate change to runoff, Water Resour. Res., 50, 9620–9629, https://doi.org/10.1002/2014WR015451, 2014.
Zhang, L., Dawes, W. R., and Walker, G. R.: Response of mean annual evapotranspiration to vegetation changes at catchment scale, Water Resour. Res., 37, 701–708, https://doi.org/10.1029/2000WR900325, 2001.
Zhang, L., Zhao, F., Brown, A., Chen, Y., Davidson, A., and Dixon, R.: Estimating Impact of Plantation Expansions on Streamflow Regime and Water Allocation, CSIRO Water for a Healthy Country, Canberra, Australia, 1–69, 2010.
Zhang, L., Zhao, F., Chen, Y., and Dixon, R. N. M.: Estimating effects of plantation expansion and climate variability on streamflow for catchments in Australia, Water Resour. Res., 47, W12539, https://doi.org/10.1029/2011WR010711, 2011.
Zhang, S., Yang, H., Yang, D., and Jayawardena, A. W.: Quantifying the effect of vegetation change on the regional water balance within the Budyko framework, Geophys. Res. Lett., 43, 1140–1148, https://doi.org/10.1002/2015GL066952, 2016.
Zhao, Q. L., Liu, X. L., Ditmar, P., Siemes, C., Revtova, E., HashemiFarahani, H., and Klees, R.: Water storage variations of the Yangtze, Yellow, and Zhujiang river basins derived from the DEOS Mass Transport (DMT1) model, Sci. ChinaEarth Sci., 54, 667–677, https://doi.org/10.1007/s1143001040967, 2011.
Zheng, H., Zhang, L., Zhu, R., Liu, C., Sato, Y., and Fukushima, Y.: Responses of streamflow to climate and land surface change in the headwaters of the Yellow River Basin, Water Resour. Res., 45, W00A19, https://doi.org/10.1029/2007WR006665, 2009.
Zhou, S., Yu, B., Zhang, L., Huang, Y., Pan, M., and Wang, G.: A new method to partition climate and catchment effect on the mean annual runoff based on the Budyko complementary relationship, Water Resour. Res., 52, 7163–7177, https://doi.org/10.1002/2016WR019046, 2016.
 Abstract
 Introduction
 Methodology
 Results
 Discussion
 Conclusions
 Appendix A: Mathematical proof of $\mathrm{\Delta}z={\int}_{L}{f}_{x}(x,y)\mathrm{d}x+{\int}_{L}{f}_{y}(x,y)\mathrm{d}y$
 Appendix B: Mathematical proof that the sum of ${\int}_{L}{f}_{x}(x,y)\mathrm{d}x$ and ${\int}_{L}{f}_{y}(x,y)\mathrm{d}y$ is path independent
 Appendix C: A fictitious example to show how the LI method works
 Appendix D: Mathematical proof of the pathaveraged sensitivity
 Appendix E: Pathaveraged sensitivity in onedimensional cases
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Methodology
 Results
 Discussion
 Conclusions
 Appendix A: Mathematical proof of $\mathrm{\Delta}z={\int}_{L}{f}_{x}(x,y)\mathrm{d}x+{\int}_{L}{f}_{y}(x,y)\mathrm{d}y$
 Appendix B: Mathematical proof that the sum of ${\int}_{L}{f}_{x}(x,y)\mathrm{d}x$ and ${\int}_{L}{f}_{y}(x,y)\mathrm{d}y$ is path independent
 Appendix C: A fictitious example to show how the LI method works
 Appendix D: Mathematical proof of the pathaveraged sensitivity
 Appendix E: Pathaveraged sensitivity in onedimensional cases
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement