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₀) as proxies for water and energy availabilities, respectively, the Budyko framework relates evapotranspiration losses to the aridity index defined as the ratio of E₀ 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=P-E$. 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₀, and n are given as

2.2 Theory of the line-integral-based method

We start by considering an example of a two-variable 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₀, y₀) to the terminal state (x_N, y_N), and z co-varies from z₀ 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$.

Figure 3A schematic plot illustrating the LI method.

Download

As shown in Fig. 3, points M₁(x₁,y₁), …, $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 total-differential 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 three-variable 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₀, 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 three-dimensional 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

1. 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₁, E₀₁, n₁) to (P₂, E₀₂, n₂) 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 integer-valued. Δ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}$$

2. 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 steady-state 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 7-year 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.

3. 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₀, E₀₀, n₀), …, (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₁, L₂, …, 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).

Table 1Summary of the long-term hydrometeorological characteristics of the selected catchments^a.

^a R, P, and E₀ represent the mean annual runoff, precipitation, and potential evaporation, all in 10⁻³m yr⁻¹. 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 total-differential 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.

Download Print Version | Download XLSX

2.4 Total-differential, decomposition, and complementary methods

To evaluate the LI method, I compared it with the existing methods, including the decomposition method, the total-differential method, and the complementary method. The total-differential 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₀, and n, respectively. Within the total-differential 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₂, P₂, and E₂ represent the mean annual runoff, precipitation, and evapotranspiration of the evaluation period, respectively; R₂^′ and E₂^′ 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₂ and E₂^′ 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 change-point 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 total-differential, 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.

Figure 4Comparisons between the LI method and the decomposition method. (a) Comparison of the estimated contributions to the runoff changes from the catchment changes (ΔR_n). (b) The decomposition method is equivalent to the LI method that assumes a sudden change in catchment properties following climate change. In this case, the integral path of the LI method can be considered as the path ABC in Fig. 3 (x represents climate factors, and y catchment properties, i.e., n) and $\mathrm{\Delta }{R}_n={\int }_{\mathrm{AB}+\mathrm{BC}}\frac{\partial R}{\partial n}\mathrm{d}n={\int }_{\mathrm{AB}}\frac{\partial R}{\partial n}\mathrm{d}n+{\int }_{\mathrm{BC}}\frac{\partial R}{\partial n}\mathrm{d}n=\mathrm{0}+{\int }_{\mathrm{BC}}\frac{\partial R}{\partial n}\mathrm{d}n={\int }_{n_{\mathrm{1}}}^{n_{\mathrm{2}}}{f}_n(P_{\mathrm{2}},E_{\mathrm{02}},n)\mathrm{d}n$, where the subscript “1” denotes the reference period and “2” denotes the final subperiod of the evaluation period.

Download

Figure 5Comparisons of the estimated contribution to runoff from the changes in (a) precipitation (ΔR_P), (b) potential evapotranspiration ($\mathrm{\Delta }{R}_{E_{\mathrm{0}}})$, and (c) catchment properties (ΔR_n) between the LI method and the total-differential method.

Download

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 semi-humid 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⁶ m². 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⁻¹, and between −35 and $\mathrm{41}\times {\mathrm{10}}^{-\mathrm{3}}$ m yr⁻¹ 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⁻¹. 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 2Comparisons of R (mm yr⁻¹), P (mm yr⁻¹), E₀ (mm yr⁻¹), and n (dimensionless) between the reference and the evaluation periods^*.

^* The subscript “1” denotes the reference period, and ”2” denotes the evaluation period. $\mathrm{\Delta }X=X\mathrm{2}-X\mathrm{1}$ (X as a substitute for R, P, E₀, and n).

Download Print Version | Download XLSX

Table 3Effects of precipitation (ΔR_P, 10⁻³ m yr⁻¹), potential evapotranspiration ($\mathrm{\Delta }{R}_{E_{\mathrm{0}}}$, 10⁻³ m yr⁻¹), and catchment changes (ΔR_n, 10⁻³ m yr⁻¹) on the mean annual runoff determined from the four evaluated methods.

Download Print Version | Download XLSX