Predicting Streamflow Elasticity Based on Percolation Theory and Ecological Optimality

How much terrestrial precipitation is used by vegetation and how much runs off, represents central issues in hydrologic science, ecology, climate change, and even geopolitics. We present a theory for the water balance to predict the fractional change in streamflow due to given fractional changes in temperature and precipitation. The theory involves a single parameter whose value is derived under the conditions of neither energy‐ nor water‐limitations and, therefore, is not an adjustable parameter. By comparison with extensive data for precipitation elasticity ϵp at global scale, we find that the theory captures the key trends of the variations of the median value of ϵp with the aridity index AI. In contrast to a shortcoming of the classical Budyko phenomenology, namely, convergence to ϵp = 4 for large AI, our theory yields a value of 2 for the median value of ϵp for all AI > 1, in accord with the data for major river basins, as well as with the median value of summaries of global and continental data sets. Incorporating in the theory the effects of annual changes in water storage leads to the ability to predict the range of observed values of the elasticity as a function of the aridity index, or its inverse, the humidity index, as well as the run‐off ratio. When changes in storage are neglected, the theory yields more accurate predictions for major river drainages than for small watersheds, particularly if the large basin spans various climate regimes and, as such, an integration over climates tends to reduce relative changes in the storage.


= + + Δ
(1) Each of the quantities in Equation 1 is typically represented as water depths accumulated or lost per year. For Q, a depth lost per year is obtained, for example, from the flux of water exiting a catchment in a river through division by the area of the drainage basin. Although the representation of water balance by Equation 1 may appear trivial, predicting E and its variability in terms of the principal input climatological variables, P and potential evapotranspiration, E p , as well as catchment storage characteristics, S, involves great complexity of the process and couplings, potentially across various length scales. In that sense, solving Equation 1 has been considered the central problem of the hydrologic sciences for about a century (Horton, 1931, as reported in Dooge, 1987NAS, 1991;NSF, 2020), although its central importance to water resources was realized well before Horton's address (1931) to AGU, having been addressed early in the 20th Century by Schreiber (1904) and Oldekop (1911). The solution will connect at a nexus of critical zone processes, thereby linking hydrology, soil science, and ecology (Eagleson, 1982;Hunt, Sahimi, Faybishenko, et al., 2023;Hunt et al., 2021aHunt et al., , 2021bNAS, 1991;Nijzink & Schymanski, 2022;Rodríguez-Iturbe & Porporato, 2007), establishing its significance in ecology, geomorphology, climate, and the history and future of life on Earth (Bonan & Doney, 2018;Nordt & Driese, 2013;Milly et al., 2005;Oki & Kanae, 2006;Vörösmarty et al., 2010;Zhang et al., 2022). One measure of the significance of a solution to Equation 1 is, for example, the result that the strongest single predictor of the net primary productivity (NPP) of plant ecosystems is the value of E (Budyko, 1974;Rosenzweig, 1968).
Solution of Equation 1 is simplified significantly when storage changes, ΔS, can be neglected, but its omission is justified at best over multi-decadal time scales (see, e.g., Gentine et al., 2012), and not at the annual scales over which data are typically reported. It has been postulated (Budyko, 1958) that the solution of Equation 1 for E(P) with ΔS = 0 can be represented as, E/P = f(E p /P), where f is a function that must approach zero in the limit E/P → 0, and unity in the opposite limit, P/E p → 0. A number of phenomenological models based on the Budyko hypothesis have been proposed in order to describe the water balance. Analytical results include those of Pike (1964), Fu (1981), Milly (1993), Choudhury (1999), Zhang et al. (2004), Yang et al. (2008), and Zhang et al. (2022). Alternatively, for a specific set of conditions, the water balance may be solved using distributed process models (Arnold et al., 1998;Henley et al., 2011;Lan et al., 2018;Martinez & Gupta, 2010;Westra et al., 2014), although their nature is not to develop an analytical solution to the water balance.
The analytical models include non-parametric ones (see, e.g., Budyko, 1958;Pike, 1964;Turc, 1954), which can only describe trends in the data, or they could be parametrized models that address as well the variability in, for example, E/P for a given E p /P. No consensus has, however, emerged regarding the controlling factors of such model parameters (Greve et al., 2015;Mianabadi et al., 2020;Zhou et al., 2015), suggesting that parameters of ad hoc formulas may not facilitate direct physical interpretation (Hunt, Sahimi, Faybishenko, et al., 2023;Reaver et al., 2022).
Distributed models also have their difficulties. Due to the complexity of processes and spatial variations among global catchments, such models generally require many parameters and wide ranges of observational data for calibration, restricting their predictive power for making accurate predictions at larger scales (Arnold et al., 1998;Beven, 1993;Moretti & Montanari, 2007). In any case, an accurate, physically based, general solution that links critical zone processes explicitly to a functional form f(E p /P) is lacking. Historically, a third approach to the water balance has been proposed to relate to ecological (or other) optimization principles (e.g., Eagleson, 1982;NAS, 1991;Rodriguez-Iturbe and Porporate, 1997). Recently, two groups, namely, Hunt, Faybishenko, Ghanbarian, et al. (2020), Hunt et al. (2021a), Hunt, Sahimi, Faybishenko, et al. (2023), , and Nijzink and Schymanski (2022) have shown that it is possible to predict essentially the same Budyko curve [the Choudhury (1999) result with parameter value 1.49] based on the principle of ecological optimality. In this paper, this approach will be adopted and applied to a direct calculation of the variability of streamflow with changing climatological variables, known as the streamflow elasticity (Němec & Schaake, 1982).
For any particular form of the function f, the associated relative changes in streamflow ΔQ/Q can be calculated as (Dooge, 1992) where ϵ p = (P/Q)(∂Q/∂P) and (1 − ϵ p ) = (P/E p )(∂E p /∂P), as long as catchment conditions do not change too severely and the climate forcing factors ΔP/P and ΔE p /E p are not too large. The quantity ϵ p is the aformentioned elasticity of streamflow with respect to the precipitation (Němec & Schaake, 1982). A valid solution to Equation 1 will provide a valid solution of Equation 2, though it may not be clear that such a solution can describe the actual behavior of catchments (Reaver et al., 2022).
In order to calculate ϵ p , we apply a recently developed theoretical approach to the water balance (Hunt et al., 2021b;Hunt, Sahimi, Faybishenko, et al., 2023; in which the effects of ΔS are incorporated explicitly in accord with Equation 1, but without a specific process model for the storage (changes). In Hunt, Sahimi, Faybishenko, et al. (2023) and , it was shown that the fundamental results previously obtained (Hunt et al., 2021b) provide a nearly perfect match for the bestfit Choudhury function for the FLUXNET data, reported by Williams et al. (2012). Using this solution, we then seek guidance for estimating the magnitude of storage variability in a given catchment by hypothesizing that our analytical formulation produces the appropriate zero storage change solution to the water balance and precipitation elasticity simultaneously. By direct comparison of the predictions for the variable storage water balance with the annual MOPEX data (that consist of 18,000 data pairs) (Gentine et al., 2012), it is then possible to estimate magnitudes of typical annual variability in storage driving the variability in E(E p , P). The magnitudes of storage change can then be used in the predictions of streamflow elasticity. As a principal associated test, we evaluate the applicability of our approach by making comparisons of the variability of elasticity ϵ p with a number of studies of continental to global extent, for which ϵ p is variously represented as a function of aridity index, E p /P, its inverse, the humidity index, or of run-off coefficient, Q/P, as well as its inverse.
Our focus on the magnitude of ΔS is shown to be consistent with the conclusion that any reduction (increase) in storage will decrease (increase) streamflow elasticity. While the approach is self-consistent, as well as generally consistent with other ideas in the current literature, we acknowledge important distinctions, as well. For example, Anderson et al. (2023) address expected effects of catchment storage properties on stream flow elasticity as follows: "low flows […] in natural rivers [depend on] inflow from catchment storage sources, such as groundwater, lakes, or wetlands (Smakhtin, 2001)," and tend to be associated with low elasticity. "High streamflow magnitudes are controlled […] by precipitation events and antecedent soil moisture conditions (Ivancic & Shaw, 2015;Slater & Villarini, 2016)," increasing streamflow response to precipitation. The authors point out that snow, with its tendency to reduce the magnitude of streamflow response, is a second example of the buffering effects of water storage. Thus, in their view, storage in different reservoirs can have widely differing effects on streamflow response.

General Theory
The primary goal in water balance research is to predict the most likely values of Q and E, as well as their distributions, as fractions of P in terms of the ratio E p /P, a formulation known as Budyko theory (Budyko, 1958). In this vein, the functions should be expressed in terms of parameters with a well-defined relationship to catchment characteristics, and not as fitting parameters to be correlated with vegetation type, drainage basin area, precipitation seasonality, etc. In principle, one would like to be able to predict not only the water balance in a given climate under current conditions, but also how the water balance will vary with changes in climate. Since the two principal climate forcing variables are E p and P, one would like to be able to predict their individual effects on stream flow, as long as climate change is not too severe. Dooge (1992) showed that the Budyko-like formulations must always satisfy Equation 2, although it does not provide information on or insight into the value of ϵ p . If we allow P and E p to vary independently, one obtains Dooge's equation, Equation 2. Thus, for example, if ϵ p = 2, then, if the precipitation increases by, say 5% and the radiant energy by 2%, climate change will contribute directly to a change in the run-off water by 8%. If vegetation characteristics change due to climate change, or if human activity changes the landscape directly, such factors may induce further variations in the streamflow. It is, of course, possible to use an equivalent formulation of the problem in terms of the thermodynamics of the phenomenon, together with Equation 1 and observations of the elasticity to predict the actual form of the water balance through integration, in which case our approach would represent a statistical mechanical formulation (Sposito, 2017).

The Percolation Approach
Plants grow in soil, a heterogeneous medium in which the local connectivity of the pores, as well as their flow conductances, which are related to their sizes, vary spatially. If the heterogeneity is strong, the morphology of 4 of 20 soil, and in particular its pore connectivity, is described by percolation theory (Sahimi, 2023). In particular, percolation theory predicts that the pore network of soil at or near its connectivity or percolation threshold, that is, the critical porosity at which a sample-spanning cluster of connected pore-space is formed, is a self-similar fractal cluster with a fractal dimension D f whose values are about 1.9 and 2.5 in two-and three-dimensional (2D and 3D) systems (Stauffer & Aharony, 1994). Such a network consists of two parts, the multiply-connected pores through which fluid flow and transport occur and is referred to as the backbone, and the dead-end or singly-connected pores. The backbone is also a self-similar fractal object with a fractal dimension D b whose values in 2D and 3D are, respectively, 1.64 and 1.87 (Sheppard et al., 1999). More importantly, percolation theory is able to accurately describe the effect of heterogeneity and pore connectedness on fluid flow and transport in porous media, and in particular transport of a solute (Ghanbarian-Alavijeh et al., 2012;Sahimi, 2012;Sahimi, Heiba, et al., 1986;Sahimi, Hugnesh, et al., 1986) in soil, and that fluid flow and transport properties, such as the permeability and the dispersion coefficients of a solute follow universal power laws near the percolation threshold (Hunt & Sahimi, 2017).
Plant roots are known to expand mostly horizontally near the surface, but are limited in the vertical extent to approximately the soil depth (see, e.g., Lynch, 1995). Since the two hydrological fluxes E and Q add up to a constant, namely, the precipitation P, the water taken by the plants is not available for making soil, and vice-versa. The percolation approach expresses the NPP, N pp , as the product of the horizontal root spreading and the vertical depth of the soil, optimizing it in terms of the hydrological fluxes, and has been able to predict the correct dependence of plant productivity on E (Hunt, 2017;Hunt, Faybishenko, & Powell, 2020;Hunt et al., 2022), as well as the dependence of soil depth on Q (Yu & Hunt, 2017a, 2017bYu et al., 2017Yu et al., , 2019. Using percolation theory, soil depth and root radial extent were derived from spatio-temporal scaling relationships by relating them to the prediction of percolation theory for the spatial and temporal scaling properties of solute transport near the percolation threshold, and the tortuosity of the optimal paths that are exploited by plant roots in order to find soil nutrients and water sources, as well as for minimizing the energy associated with media deformation. Solute transport is highly important to soil development, since its role in the fundamental process of chemical weathering that underlies soil formation is critical (Egli et al., 2018;Yu & Hunt, 2017a). In the case of vegetation growth, the percolation scaling based on soil pore network imposes limitations on plant root growth. The fundamental length scale is the median particle size, which represents the distance between neighboring pores, while the ratio of that length scale and the related hydrologic fluxtranspiration for plant growth and subsurface run-off for soil development-defines the dominant time scale. Based on such considerations, N pp is related to the morphology of soil through the following relation (Hunt et al., 2021b), Here, D f = 91/48 ≈ 1.9 is the aforementioned fractal dimension of 2D percolation networks, proposed (Hunt, Faybishenko, & Powell, 2020;Hunt et al., 2022) to describe the plants' root network, and 1/(D b − 1) ≈ 1.15, if we use the fractal dimension of the 3D percolation backbone. The reason for the distinction between the fractal dimension D f of 2D networks in one case and the backbone of 3D networks, represented by D b , in another case lies in the 3D nature of the soil network, as experienced by a solute that is being transported through the soil, as opposed to the 2D nature of the soil-root network, since the root tip extension has been confined to be within the upper boundary (Gentine et al., 2012;Lynch, 1995) and, therefore, it is essentially a 2D structure. Optimizing N pp as given by Equation 3 with respect to E, that is, setting ∂N pp /∂E = 0 and solving for E, yields where d is the spatial dimension, as both D f and D b dependent on d. If we use D f (d = 2) ≈ 1.9 and D b (d = 3) ≈ 1.87, we obtain E = 0.623P. It is important to emphasize the two approximations in Equation 2, namely, substituting evapotranspiration E for transpiration, and the subsurface component Q s of the run-off for the run-off Q. The approximation is rather accurate since transpiration is typically about 70% of E (Schlesinger & Jasechko, 2014), while the subsurface component of run-off is a similar fraction of the total run-off (Jasechko, 2019). Invoking such approximations does, however, admit an important input into the variability between watersheds of the actual streamflow, which is assumed to be equal to Q, which we do not specifically include.
The question of water limitations is addressed by assuming that the optimization is valid only for the vegetated portion of the landscape, and that the fraction is proportional to P/E, an assumption that has been verified by comparing with data (Yang et al., 2009), while we assume that evaporation in between the plants is 100%. The question of energy-limitations is addressed by assuming that for P > E p (potential evapotranspiration) the optimization could only apply to that P = E p and the remaining P runs off. The prediction for the water balance has been compared with both long-term and annual data (Hunt, Faybishenko, Ghanbarian, et al., 2020;Hunt et al., 2021b).
While we can make specific predictions for water balance, we must do so for each range of A I values separately, since we derive distinct results for the separate ranges. The water balance functions derived in this way are given by Note that the numerical values of α(d) and their dependence on the spatial dimension d, α = 0.623 and 0.813 for, respectively, d = 2 and 3, are due to whether the plants' root system is mainly confined to the very shallow subsurface, that is, a 2D system, or is more nearly isotropic in the case of isolated plants, competing for water in 3D. Note also that the value of 0.623, obtained by assuming a very shallow soil and root system and no limitations on energy or water, yields very nearly the global average E/P = 0.632, as given by the mean value of the data in eight studies that reported E/P to be ranging from 0.61 to 0.67 (Hunt, Faybishenko, Ghanbarian, et al., 2020).
It is straightforward to verify that the above two expressions both yield the constant 1 − α(d) = 0.377 at A I = 1. This constant is of considerable importance for interpreting the dependence of ϵ p on the run-off coefficient R, since the functional dependence will change at this value of R, not as obvious a watershed value as A I = 1 is for climate representations of the elasticity. The slope of the functions-their derivative with respect to E p -is, in general, not continuous, implying a profound difference between the water-limited (A I > 1) and energy-limited (A I < 1) regimes, which is in agreement with the qualitative differences in the response of watersheds to climate change in the two regimes (see below).
Next, we calculate the precipitation and energy elasticities, defined above, for our particular Budyko-like function. We first use the expression for Q corresponding to A I > 1 given in which, given the expression for Q, is equivalent to, This implies that in Dooge's formulation, we obtain a logarithmic derivative of Q with respect to P, and ϵ p = 2, that is, a precipitation elasticity of 2 and energy elasticity of −1. We can also use the expression for Q for A I < 1 to obtain which also satisfies Dooge's complementarity, but not with simple coefficients that are independent of A I .
We also note that since α is explicitly dependent upon vegetation strategy, it can be used to represent the effects of vegetation changes on the water balance. Since a more isotropic root system is likely to be deeper, an increase in α would also be consistent with an increase in the rooting depth, or active soil layer. Such a procedure would add the terms The parameter α is sensitive to rooting depth, such that increasing it tends to increase the value of α. If, for A I > 1, increasing aridity due to either a decrease in P or an increase in E p should cause rooting depth to increase, that would result in a vegetation contribution that enhances the contribution due to either P or E p alone and, hence, a value ϵ p > 2. Consider, for example, a situation in which P increases in a dry climate, adding grasses to a landscape that had none. In such a case, it would be conceivable to reduce α(d) from 0.813, its value for d = 3, to 0.623, which is for 2D media, with the result being an increase in ΔQ/Q (about equal to) dQ/Q of 0.187/0.19 ≈ 1, in addition to the direct contribution of a factor of 2. Conversely, the loss of grass due to overgrazing and its replacement by woody shrubs with very deep roots could cut streamflow in half, that is, −0.187/.377 ≈ −0.5, even in the absence of climate change.
In the remainder of this paper, we consider only the 2D value of α(d = 2) = 0.623, partly because of its close correspondence with the global mean value of E/P, and partly because the agreement with experiment on the streamflow elasticity is not enhanced by considering potential variability in α(d) as well. Moreover, the predicted median value of ϵ p is the same for A I > 1 for either choice of α(D). Although uncertainty regarding the distinction between evapotranspiration of forests versus grasslands exists, woody plant encroachment (WPE) in regions historically associated with grasslands has often been cited as a reason for diminished streamflow. Further, many factors may influence the actual change in streamflow (Huxman et al., 2005), yet the mean of the run-off coefficients from 9 studies of WPE cited in their Table 1 is 70% higher with "low woody plant" concentration than with "high woody plant" concentration, rather close to our estimate of 100%, if the change of root fractal dimensionality is from 3D to 2D in all plants in the region.
We point out that it would be difficult for the present analysis to account for a separate reduction in vegetation coverage due to increasing aridity, since an increase in A I is already assumed to decrease vegetation coverage, represented as 1/A I . Evidence reported by Yang et al. (2009) suggests that such an onset of vegetation loss occurs at A I > 1. Nevertheless, although an increase in A I at a given location for A I > 1 does not produce an explicit contribution to the precipitation elasticity, it would seem to be consistent with an additional contribution to stream flow decrease.
In order to understand better the data for the annual streamflow precipitation elasticity, we address the effect of the potential change ΔS in storage S in the 2D systems. Writing and differentiating it yields, which means that (where 2.65 ≈ 1/0.377) which does not follow Dooge's theorem, but approaches the equality in the limit ΔS → 0. The final term, a storage elasticity, will be neglected in the following, but was included for completeness.
Adding the storage term to Equation 4 and repeating the derivation that led to Equation 8 yields, Since, we assume that for A I < 1 energy limits the productivity, whereas water does so for A I > 1, we make an analogous argument for storage. Thus, the maximum |ΔS| for energy-limited systems is a small fraction of E p , while for water-limited systems, it is the same fraction of P. For a reasonable estimate of year-to-year changes in storage, we estimate the fraction to be 10%. In the Great Lakes watershed values of ΔS of about 7% have been reported (Niu et al., 2014), while in western USA values as high as 30%-40% have been cited (Garcia & Tague, 2015). Focusing on the precipitation elasticity, and keeping in mind that the sign of ΔS may be positive or negative, we find that In Equation 15, we assumed ± 10% change in storage in terms of E p or P, as discussed above. Note that the contrasting assumptions regarding the relevance of E p and P to storage in humid and arid climates, respectively, leads nevertheless to results that predict dependence of ϵ p on A I , which is increasing for either A I < 1 or A I > 1. In the cross-over from the humid to the arid regime, however, there is a sharp drop in ϵ p .
If positive changes in storage are as likely as negative ones, our results that neglect storage should be interpreted as the median values for the precipitation elasticity. For the same absolute value of storage change, however, positive changes in storage predicts a larger increase in ϵ p than do negative changes in storage in their reduction of ϵ p . Thus, the median value of ϵ p is not the mean value, and the latter expected to be somewhat larger than the former.
Some authors use the inverse of the aridity index, P/E p , called the humidity index, as the independent variable in presentation of data for ϵ p . Some authors also use the run-off coefficient, Q/P (or its inverse) to present data for ϵ p . Note that our result for ϵ p = 1/[1 − α(d)A I ] for A I < 1 is also equal to P/Q, as is, in general, any equation for E that is a linear function of P and E p .

Data Sets
In order to compare our predictions for ϵ p with observation, we restrict our search to values obtained non-parametrically, while seeking data from the largest compilations currently available. In addition to the large sources, we present graphically results obtained from the relatively small database of Zhou et al. (2015) with values of non-parametrically obtained streamflow elasticities for mainland China. It is known that ϵ p may depend on the details of the precipitation anomaly, as well as on the time period associated with measurement . With one exception (Deusdara-Leal et al., 2022), we have considered only data sets for which the reported ϵ p are obtained over annual periods. Similarly, aside from Deusdara-Leal et al. (2022), the authors, whose data we accessed, made no distinction between positive and negative anomalies in P. The climatic data of Sankarasubramanian et al. (2001) were generated from PRISM [see Daly et al. (1994)], while the streamflow data derived from the Hydroclimatic Data Network [see Slack et al. (1993)]. Because data for only 5 of the 17 USA regions mentioned were explicitly given in Sankarasubramanian (2001), we also acquired the associated dissertation of Sankarasubramanian (2001). However, no additional data were found. Thus, what we could digitize included only the published figures, covering regions 1 (New England), 3 (the southeast), 10 (the Missouri River drainage), 12 (the western Gulf of Mexico), and 17 (the Pacific Northwest), with median (mean) ϵ p = 1.73 (1.82), and (excluding a single data point of A I > 9) a maximum A I of 4.5. The entire data set included 1291 catchments with at least 20 years of records and areas exceeding 129 km 2 . For each catchment, the reported ϵ p , independent of precipitation, was the median of annual measurements, that is, where, for example, ̄ = mean value of P. In a separate work, Sankarasubramanian and Vogel (2003) presented data for ϵ p as a function of inverse run-off ratio P/Q, which we also digitized. The data are, nevertheless, from the same source. The data of Harman et al. (2011) included 405 catchments across the USA obtained from MOPEX [see Duan et al. (2006)]. Catchments of sizes 80 km 2 to more than 10,000 km 2 , for which streamflow data were available for at least 25 years, were used. Our digitized data from Figure 1 (top panel for total run-off elasticity) of this source had a median elasticity of 2 and a mean elasticity of 2.12 (the authors reported variously a mean of 2.10 and 2.08). The data reflected a maximum A I of 4.2, although the authors reported an absolute maximum A I of 8. Harman et al. (2011) did not use a non-parametric determination of ϵ p , opting for a procedure employing the Ponce and Shetty (1995) water flux and balance model. They also split the run-off into fast and slow components, but reported simultaneously a combined elasticity, which is what we used.
We digitized the figure on page 12 of Chiew et al. (2006) for ϵ p as a function of Q/P with sites on every continent except Antarctica, as well as on islands, such as Ceylon, New Zealand, and Iceland. The digitized data had a median ϵ p = 1.68 and a mean ϵ p = 1.74. The criteria for selection was possession of concurrent sources for monthly streamflow (Peel et al., 2000), precipitation, and temperature (Global Historical Climatology Network, GHCN). This stipulation allowed calculation of the elasticity non-parametrically for 532 catchments with time periods from 23 to 64 years and areas ranging from 100 km 2 to 76,000 km 2 . For each catchment, the median of the annual elasticity values was determined according to the procedure of Sankarasubramanian et al. (2001) described above.
The digitized "unimpaired," that is, dam-free, Australian data with source Chiew (2006) for ϵ p as a function of HI (their Figure 8c, panel 2) and as a function of Q/P (their Figure 8a panel 2) had a median value of 2.43 and a mean value of 2.57, as well as a maximum A I of 6 (one value), with the remaining data confined to A I < 5. The catchments ranged in size from 50 km 2 to 2,000 km 2 . Note that we avoided Figures 8a and 8c of panel 1, since the ϵ p values given there were model-dependent.
The major river basins compiled in the supplement to Tang and Lettenmaier (2012) included 25 catchment areas between 80,000 and 100,000 km 2 , 26 catchment areas that exceed 1,000,000 km 2 , and the remaining ca. 140 with areas between 100,000 km 2 and 1,000,000 km 2 , covering roughly half the world's continental area. Tang and Lettenmaier (2012) thus provided a data set complementary to all the remaining data. In this calculation/compilation of ϵ p , the 194 largest river basins were chosen [referred to as Simulated Topological Network (STN-30p); see Vörösmarty et al. (2000)]. The authors did not report the aridity index for the basins. Since Liu et al. (2018) reported A I data for 28 major river basins, we contacted one of the co-authors, Yongqiang Zhang, who provided 400 additional A I values that had been obtained using the methods described in Lan et al. (2018), from which we were able to generate 125 data pairs (ϵ p , A I ) after accounting for different names of the same rivers. Zhou et al. (2015) calculated ϵ p for basins within China by non-parametric methods, as well as by the Budyko method, but represented ϵ p as a function of A I for the non-parametric methods. According to Zhou et al. (2015), the three sites with very low values of ϵ p are cold sites. Overall, the regions that produced the smallest values of ϵ p in the very humid southeast, or cold and dry climates in the far northwest.

Alternate Model
When assessing the correspondence with actual data, it will be illustrative to apply also an existing phenomenological model for comparison, of which we choose Budyko's (1958) function. The Budyko (1958) function is venerable and still often used (see, e.g., Berghuijs et al., 2020;Mianabadi et al., 2020;Zhang et al., 2022), although it has also been subject to criticism (Berghuijs et al., 2020;Greve et al., 2015;Reaver et al., 2022), Its lack of parametrization restricts is versatility. One might also object to its origin as a synthesis of two previously existing phenomenological models (Oldekop, 1911;Schreiber, 1904) (Duan et al., 2006) for 50-year averages of E/P, showing that the Budyko function gives a good accounting for (E/P) as a function of (E p /P) throughout the range of climates. Gentine et al. (2012) argued that any discrepancy between the Budyko curve and the long-term average data could be attributed to experimental error of about 10%. Moreover, the purpose of Zhang et al. (2022) was to parametrize the Budyko function in order to deal with variable storage changes.
In the present context we employ the Budyko function only as a benchmark prediction for comparison with our own predictions for precipitation elasticity of Q, as Zhang et al.'s (2022) inferred changes in E/P are quite different from what we will report here. Ultimately, it will be the comparison with direct field data which should guide one's choice of theory. Sankarasubramanian et al. (2001) found that median ϵ p values of individual catchments ranged from 1 to 2.5 over the entire USA with the exception of a small (cold) portion of Montana and North Dakota with lower values. These values compare well with our prediction, namely, that over the range 0 < A I < 1, the median elasticity should vary from 1 to 2.65 and for A I > 1, the median elasticity should be 2. Note that in only 0.7% of the range of A I values, from 0.964 to 1, is our predicted median elasticity greater than 2.5, while its absolute maximum is 2.65. In regard to the data for North Dakota and eastern Montana, studies of streamflow elasticity in catchments with significant snow cover often reveal ϵ p values near 1 and even smaller (e.g., Berghuijs et al., 2016;Khan et al., 2022;Sankarasubramanian et al., 2001;Zhou et al., 2015). This has been interpreted as a buffering effect of snow ; clearly, time lags between P and run-off may be considerable, while, on shorter time scales, run-off may most strongly be positively correlated with solar energy, rather than P. Xing et al. (2018) investigated the precipitation and energy elasticities of 35 drainage basins in China. They found mean values of P and E p elasticities to be 1.972 and −0.972, respectively, satisfying the Dooge (1992) relationship. The range of P elasticities in their catchments was from 1.24 to 2.45; the E p elasticity values again satisfying the Dooge (1992) relationship. The geographical distribution of their elasticities produces the smallest values in high altitudes with drought conditions, where E p is small, and in the far southeast where P is largest.

Data Summaries
For the Kuye River basin, Zheng et al. (2021) calculated pre-anthropogenic and post-anthropogenic values for the precipitation and energy elasticities, using three different phenomenologies for the water balance, rather than non-parametric methods. Only in the pre-anthropogenic period (1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966) do the values obtained using all the distinct phenomenologies agree closely; in that case, ϵ p = 1.95 ± 0.04. Other authors (e.g., Chiew et al., 2006;Zheng et al., 2009) have discussed [or, Zhou et al. (2015), given data consistent with] the tendency for the precipitation elasticity to increase with increasing aridity for A I < 1 and then, particularly the median value of its spatial variability, to tend toward independence of A I when A I > 1; our prediction for A I > 1 when ΔS = 0, is ϵ p = 2, which appears to be the most commonly observed overall value of the P elasticity. The prevalence of precipitation elasticity values near 2 is confirmed in a study of 521 catchments worldwide by Chiew et al. (2006). They presented a table that we reproduce here (see Table 1). Note that the range of ϵ p values is defined by Chiew et al. (2006) as lying between the 10th and the 90th percentile, while the mean of the ϵ p values at these two percentiles is greater than or equal to the median in every climate category except 1.

Inferences Regarding Storage Effects
A number of authors (see e.g., Berghuijs et al., 2016, and references therein) have noted that annual changes in storage should have significant effects on the water balance, as well as on precipitation elasticity of streamflow. Drawdown of storage reflects, in part, extra water available for streamflow, which can reduce the effects of negative fluctuations in precipitation over the same time period. While an increase in storage reduces immediate run-off, it can also make streams more sensitive to precipitation fluctuations over longer times. Since our investigation uses data from other sources that did not provide water storage changes, we are unable to use the storage variable predictively, but only diagnostically. Thus, we sought guidance from two sources: (a) references to storage changes as a function of precipitation for specific catchments with known climates, and (b) a data summary for the water balance that could allow us to estimate storage variability from our theoretical predictions. For the former, we refer to Niu et al. (2014) who found in Great Lakes watersheds with A I near 1 annual variability of ΔS ca. 7%P, while Garcia and Tague (2015) give specific estimates of storage change in the western USA From Oregon and California to Colorado, annual variability of ΔS could range from about 20% to 30%-40%P. For the latter estimates, comparison of our water balance prediction with the MOPEX annual water balance data published by Gentine et al. (2012) (After Chiew et al., 2006) balance extracted from MOPEX (Duan et al., 2006) by Gentine et al. (2012), which, in Figure 1 are compared with E/P, as obtained from 1 − Q/P from Equation 4 above and invoking the effects of storage through Equation 1.
As can be seen in Figure 1, storage changes ΔS with values −30% < ΔS < 10% (blue and red solid curves), as applied to a median water balance obtained from the optimization of NPP for root systems embedded in 2D, appears to be sufficient to constrain the water balance data over the entire range of A I values. In fact, for A I > 1, this range of storage changes overestimates the water balance variability both above and below the ϵ p = 2 solution. Thus, we propose that the same limits on ΔS should, in most cases (particularly within the conterminous USA), serve to constrain the elasticity data. As discussed in our theoretical formulation, this percentage is referenced to E p (P) for A I < 1 (A I > 1).
For large A I , it is likely that storage changes proposed to be as large as (−30%) are unrealistically high, since total shallow subsurface storage may be a small fraction of P. Thus, we have also shown a cut-off at E/P = 1 when the predicted E/P exceeds 1. The data are in general accord with this suggestion, being confined to the narrower range of values than predicted in the limit of large A I , both above and below the assumed median.

Comparison With Streamflow Elasticity Data
In Figure 2, we compare the values of ϵ p reported by Tang and Lettenmaier (2012) for the world's major river basins and Zhou et al. (2015) for selected rivers in China, with our predicted dependence, but using storage changes of only +8% and −10% in Equation 14. The results for ϵ p for major river basins (typical sizes in the hundreds of thousands of square kilometers) fall comfortably within the predicted limits, but the data from China exhibit a noticeably wider scatter than is predicted by these storage changes. While a decrease in relative storage magnitudes with increasing basin size is expected from general arguments regarding the summation of, for example, uncorrelated random variables (at length scales longer than typical atmospheric Rossby wavelengths, i.e., 500 km, one often finds negative correlations), it may also be a consequence of the different means Tang and Lettenmaier (2012) used to generate ϵ p . Clearly, the phenomenology of Budyko (1958) is not well adapted to the data shown, although incorporation of ΔS variability would allow a more realistic appearance (yet without bringing the predicted median value into closer correspondence). The prediction of a median value of ϵ p = 2 for A I > 1 Figure 1. A demonstration that the 2D prediction of percolation theory together with annual storage changes of 10% P (E p , for A I < 1) and −30% P (E p for A I < 1) bounds the MOPEX annual water balance data accurately for A I < 1. For A I > 1, the data scarcely approach these bounds; moreover, −30% P exceeds significantly the usual Budyko limit imposed for longer time frames. MOPEX data from Gentine et al. (2012) were digitized by Hunt et al. (2021b) and re-used here. Note the very close correspondence between our optimal carbon assimilation predictions (Hunt et al., 2021b) and the optimal carbon profit model of Nijzink and Schymanski (2022), both of which predict a water balance that is very close to the mathematical model of Choudhury with n = 1.49, although Nijzink and Schymanski (2022) include also a "low optimum" that corresponds to n = 1 as well as our result with a 10% storage increase.
looks to be in accord with the Tang and Lettenmaier data, and is in quantitative agreement with the conclusions of Chiew et al. (2006) mentioned above.
In Figure 3, we present a comparison of predictions of Equation 14 for ϵ p with the corresponding data in three regions presented by Sankarasubramanian et al. (2001). The three regions, 1, 3, and 12, correspond to New England, the southeast, and western Gulf of Mexico, respectively, which document the kind of response expected with increasing A I . Note the increase in the median value of ϵ p with increasing aridity index for A I < 1, but no apparent increase for A I > 1 (where the data are noticeably more symmetrically distributed about ϵ p = 2), in  (2012) and Zhou et al. (2015) and using the storage changes mentioned. The "Predicted elasticity 2D" curve refers to zero storage change, interpreted in the text to be a median value. For reference to traditional theoretical descriptions, we include the phenomenology of Budyko (1958) as well. agreement with our predictions. Thus, what appears at first glance to be purely an artifact of data collection provides additional evidence for distinct dependence of elasticity on A I for A I > 1 and A I < 1. Figure 4 presents a comparison of our predictions for ϵ p with data for the Missouri River basin and the Pacific Northwest. The former data set is concentrated mainly in the region A I > 1, but, owing to the strong rainshadow effects of the Cascade mountain range, the latter data set spreads across both humid and arid regions. Again, the predicted tendency for the median value of ϵ p to stay constant at A I > 1, but to increase with increasing A I for A I < 1 seems to be borne out, although here, the cross-over may occur at a lower A I (higher humidity index) than predicted.
In Figure 5, we present the MOPEX data of Harmon et al. (2011) and the Australian data of Chiew (2006) for ϵ p as a function of humidity index H I = 1/A I . Equation 14, with storage changes ΔS = ±13%, provides reasonable upper and lower bounds for the observed data, significantly smaller than the range of storage changes chosen to constrain the original MOPEX water balance data reported by Gentine et al. (2012). While there is no additional confirmation of a peak in ϵ p at A I = 1, the predicted changes in slope and curvature appear to be verified. The median ϵ p value at large A I appears to be larger than 2, although the median ϵ p value for one of the (entire) data sets is reported to be slightly less than 2 (Harmon et al., 2011). Chiew et al. (2006) do report a median elasticity larger than 2 for Australia, however-the only continent with such a feature. Although not explicitly shown here, the Budyko (1958) prediction would perform somewhat better here. A later figure presents the four main data sets together. Figure 6 inspects data for ϵ p that were reported as a function of inverse run-off ratio, P/Q, rather than as aridity index, A I . The data (Sankarasubramanian & Vogel, 2003) originate from the same study, Sankarasubramanian et al. (2001), as represented in Figures 3 and 4, however. Thus, one might expect a similar validity of prediction here. However, a somewhat wider range of storage changes needed to be applied, in order to generate a comparable correspondence with the data. A possible explanation for the apparent added variability is the sensitivity of the run-off ratio itself to variables that are unaccounted for. The vertical line in the curve appropriate for ΔS = 0 is predicted for A I = 1 to fall at P/Q = 1/0.377 or 2.65. The current representation in terms of P/Q, thus, obscures the particular relevance of A I = 1 by distributing the predicted sharp cross-over in terms of A I over different P/Q values from each distinct value of ΔS.
In Figure 7, we represent the dependence of ϵ p on the run-off ratio, Q/P, as reported by Chiew et al. (2006) from a global study (a) and an Australian study (b). Our predictions for various storage changes are given. Each branch  . As in Figure 3, the limits on ΔS implied by the comparison of the water balance prediction with the MOPEX data summarized by Gentine et al. (2012) are applied. The Budyko prediction is again applied for comparison.
of the curve is cut off either at the margin of the figure, or at A I = 5, the largest A I value reported by the authors. Predictions for each ϵ p change sharply at A I = 1, but at different values of Q/P. Several features of the data are picked up by our theoretical framework, including, (a) if the change in storage ΔS could take on any value, the predictions would ultimately trace out the entire curve of ϵ p = 1/(Q/P), our result for small A I (a result considered an upper limit by Sankarasubramanian & Vogel, 2003). However, for the given data, this bound becomes irrelevant for ΔS greater than about 10%, even though such large storage changes can still be relevant for A I > 1.   and Vogel (2003), which are, however, derived from the same data set that was published in Sankarsubramanian et al. (2001).
Here, it appeared that a somewhat larger positive ΔS value of 15% (rather than 10%) was needed to constrain the large majority of the data. Note that, as observed by Sankarasubramanian and Vogel (2003) and implied as well in the publications of Chiew (2006) and Chiew et al. (2006), we predict that all data should fall below the line given by the equality of inverse runoff ratio ( ) and ϵ p . Clearly this prediction is violated. However, the implied result that Q not exceed P is also violated, and the discrepancies with the predictions of ϵ p are largest where Q > P. In constructing these predictions, the same change in storage was entered into both the run-off ratio and the calculation of ϵ p .
(b) Cutting off the predictions at A I = 5 is in rather close correspondence to the data. (c) The sharp change in behavior of ϵ p at A I = 1 shows up in the data near ΔS = 8%. (d) For small run-off ratio, changes in storage can be important in extending the range of run-off ratio values, for which the theory applies, to near zero. It should be noted that the Budyko approximation, without addition of changes in storage, does not particularly capture the trends of the data.
In Figure 8, we compare the predictions of Equation 14 for ϵ p as a function of humidity index with all the major data sets simultaneously. We have used the bounds on ΔS, which were found to be compatible with the Figure 7. Here, the annual data for ϵ p as a function of the run-off coefficient , reported by Chiew (2006) and Chiew et al. (2006) are presented, together with the Budyko (1958) curve and our predictions for ϵ p . In this comparison it was necessary to include ΔS values as large in magnitude as (-) 55%, much larger than in any other study. Each branching curve is cut off either at the margin of the figure or at the largest A I value found in the data, which is 5, as given in the text. The Budyko (1958) prediction is given for comparison. Figure 8. Comparison of all major data sets for ϵ p as a function of Humidity Index. Bounds on the storage change are taken to be ΔS = 10% and −30%, the values compatible with the MOPEX data for the water balance predictions. The Budyko (1958) curve is shown once again for comparison. MOPEX data obtained from Gentine et al. (2012), which are +10% and −30%. Of all the data sets, only that of Sankarasubramanian et al. (2001) has a significant number of points (ca. 25), lying significantly outside the boundaries associated with the storage changes given. The zero-storage change prediction, unlike the Budyko curve, remains near the middle of the range of data points, consistent with a relatively universal result from the Chiew (2006) data summary, giving values of nearly 2 for the median elasticity in almost every climate class. The predictions are much less sensitive to changes in storage at large values of the humidity index, in agreement with the data. By following separately each of the two predicted peaks in elasticity at A I = 1 and large A I , each of the four data sets shows evidence of the relevance of the cross-over in functional form for ϵ p predicted at A I = 1. The Tang and Lettenmaier (2012) results exhibit, as expected, an overall narrower range of ϵ p , particularly at very low values of the humidity index, than all other sources. Thus, a much smaller range of relative storage change values is sufficient to account for the variations in elasticity for very large river basins than for the smaller catchments typically chosen for investigation in large studies.
Finally, consider the study of Deusdara-Leal et al. (2022) on streamflow in southeastern Brazil. They reported data for Q, A I , and the streamflow elasticity ϵ p . We first show that our 2D prediction for the water balance is accurate for their streamflow indices; the comparison is shown in Figure 9a. In Figure 9b our prediction for ϵ p (i.e., lines with specified slopes), as calculated for zero storage change, are compared with annual data for the ΔQ(ΔP). The predictions (blue and red lines) correspond to the expected streamflow precipitation sensitivity for the corresponding aridity index values, which are near 1, but neglecting effects due to storage change. When A I is very close to 1, the elasticity is very sensitive to actual A I . An obvious result from the comparison is that our predictions for the elasticity underestimate the actual values, since a significant fraction of the data fall outside our predictions (between the red and blue lines in the figure). Interestingly, however, a large share of the data comes from multiyear correlations. The precipitation was also anomalously low for almost all the regions for the entire period. Thus, cumulative effects of continued drying of the CZ may accentuate the values of the streamflow elasticity over time. Notably, however, the one river basin that follows our minimum slope calculation has undergone forest restoration, helping to stabilize soils as well as water storage.

Discussion
Use of a new theoretical study to compare streamflow elasticity data in different datasets has made it possible to develop new understanding. The effects of storage (change) should be felt in streamflow elasticity, since water going into (coming out of) storage reduces (increases) streamflow. By comparing our predictions for the water balance with and without storage changes, we were able to estimate reasonable bounds on the magnitude of interannual storage changes that could impact the streamflow elasticity, particularly in the USA, but also globally (to a lesser degree).
We found that there may be advantages to representing precipitation elasticity as function of the run-off ratio, since that allows simultaneous assessment of the relevance of storage changes and aridity index, thus a measure of the diminishing relevance of large magnitude storage changes at large A I values. This finding is in accord with the understanding that storage in arid regions is typically small. In our interpretation, a relatively universal result that the median value of ϵ p for A I > 1 is 2 is a consequence of the relevance of the case of zero storage change to the median value of the precipitation elasticity.
Our interpretation of the relatively small range of values of ϵ p for any particular value of the aridity index in the world's great river basins is that, it is due to the rarity of having the same sign of subsurface storage changes over the very large areas of such basins, particularly when they span multiple climate regimes. In a related matter, Zhang et al. (2021) found that catchment elasticity values tended to increase with time over sub-annual time scales, but at times exceeding 1 year, tended to saturate at values between 1.3 and 2.3 with a median value near 1.8, not very different from the summary of Chiew et al. (2006). Clearly, the Australian catchments studied by Chiew (2006) have a greater sensitivity to precipitation anomalies (larger values of ϵ p ) than the worldwide selection of catchments in Chiew et al. (2006). If our theoretical model is a guide, this is associated with a greater importance of positive storage changes in Australia, as compared with negative storage changes globally. This result appears to be counterintuitive, but it is in accord with the study of Deusdara-Leal et al. (2022) for southeastern Brazil. Perhaps when drought is persistent, a minimum storage is reached and storage changes are mostly positive, a negative correlation based on antecedent anomalies. We note that our failure to capture the magnitude of the southeastern Brazil elasticity values is then likely due to a neglect of potential storage changes.
We suggest that the small values of precipitation elasticity in snowy catchments is related to the phenomenon of snowmelt. This has three simultaneous implications: Absorption of solar energy through warming the snow and the phase change of melting should tend to reduce effective value of A I . It is known (Dooge, 1987) that the precipitation and temperature elasticities sum to 1. On sub-annual time scales, snowmelt, and particularly ice melt, tend to make the streamflow response to solar energy positive, with consequence that the response to precipitation must trend smaller and more nearly negative. Finally, as other authors have also noted (e.g., Sankarasubramanian et al., 2001), snow can be regarded as buffering streamflow response or simply as another form of water storage.
The data reported by Zhang et al. (2022) imply that accounting for storage results in nearly constant E (E p /P) for A I > 1, and very little difference from the original Budyko function for A I < 0.5, whereas our predictions imply a strong dependence of E (E p /P) for small A I , as well as a strong dependence of ϵ p on storage changes, except for small A I . Although not the focus of the present paper, these two perspectives are certainly very different.
It is interesting that our prediction of the water balance is highly successful for the Brazilian data, but that we underestimated the elasticity significantly. We propose without proof that the long-term drought in this region could explain the underestimation of the elasticity in our prediction. We suggest that the close correspondence between the particular results of Nijzink and Schymanski (2022) for their water balance calculation through optimization of net carbon profit with highest evapotranspiration (consistent with Choudhury's n = 1.49) and the results of our optimization of carbon assimilation implies greater reliability of each study individually.

Conclusions
We have applied the fundamental definition of precipitation elasticity of streamflow to an existing model of the water balance. The model was formulated in terms of an ecological optimum, that is, the ecosystem we find to be present at any particular location, is the one that is most successful at drawing down atmospheric CO 2 to grow and reproduce. By setting the derivative of the NPP with respect to evapotranspiration equal to zero, it was then possible to derive the actual values of the principal hydrologic fluxes, evapotranspiration E and streamflow Q, as functions of potential evapotranspiration E p and precipitation P. Extension to water-limited and energy-limited ecosystems involved restricting the procedure to the area covered by vegetation and the fraction of precipitation equal to potential evapotranspiration, respectively. Previous work documented the suitability of this procedure for describing the water balance itself. It was suggested, however, that such an abrupt cross-over in functional form at the boundary between water and energy limited environments could be unphysical and impair the success of future research based on this hypothesis.
We have shown that the abrupt boundary with its sharp change in slope of the precipitation elasticity as a function of aridity index at A I = 1 is reflected in all the major data sets that report ϵ p as a function of either aridity index A I , or its inverse, the so-called humidity index. This cross-over is harder to identify in cases where the precipitation elasticity is graphed as a function of run-off ratio, Q/P, or its inverse.
We predicted that for A I > 1 the precipitation elasticity should have a median value of 2, whereas for smaller A I the median value should rise gradually from 1 to (at A I = 1) 2.65. Storage changes will affect values of elasticity determined on annual time scales, adding significant spread to these values. Guided by a test of the relevance of specific ranges of storage values to the water balance itself, we were able to predict ranges of elasticity that were in good agreement with the experimental results. The tested relationships performed much better than the popular phenomenology due to Budyko, though it is not clear how much that traditional phenomenology might be improved by making it compatible with annual storage changes as in Zhang et al. (2022).